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The Static Stability of Half-Bubbles 


By W. J. DeBONTE 
(Manuscript received May 16, 1972) 


A variational model is used to calculate the static stability limits and 
equilibrium properties of ‘“half-bubbles,”’ magnetic domains residing on 
one surface of a magnetic bubble material platelet. Stability is achieved 
through the presence of gradients in the domain wall energy density and/or 
saturation magnetic moment. The model evidences two distinct types of 
instability behavior separated by a critical value of the wall energy density 
gradient. Unlike the standard cylindrical domain, the half-bubble has a 
minimum stable value of the ratio of domain diameter to height. 

The half-bubble 1s shown to possess a number of properties which make 
at potentially useful for device applications. It is self-biased by its closure 
wall and in some cases is stable in zero external bias. Bias margins are of 
the same order as those for the standard cylindrical domain. It stabilizes 
on only one platelet surface, and its properties are independent of both 
material thickness and of minor irregularities on the second surface. In 
addition, its structure may be advantageous in avoiding the undesirable 
properties of hard bubbles. 


I. INTRODUCTION 


The current interest in cylindrical magnetic domains and their device 
applications has stimulated the search for such domains in a wide 
variety of materials produced by a number of growth techniques. In 
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some samples, Bobeck’ has observed domains of reversed magnetization 
which apparently do not penetrate through the entire thickness of the 
sample (Fig. 1). These domains, referred to as ‘“‘half-bubbles,” have 
four features which make them potentially attractive for device applica- 
tions. First, they contact only one surface of the platelet in which they 
reside so that their properties should not be as critically dependent on 
surface preparation as those of the usual cylindrical domains. Second, 
as shown in this paper, their properties are independent of platelet 
thickness. Third, this article also shows that their top closure wall 
produces a self-biasing effect which makes external biasing unnecessary 
in some cases. Finally, Bobeck et al.” have observed in two-layer films 
that domains having a closure wall do not exhibit the undesirable 
properties of hard bubbles, and it is a reasonable presumption that this 
property also pertains to half-bubbles. 

Bobeck* has suggested that half-bubbles may be stabilized by the 
presence of gradients normal to the material surface in one or more 
of the material characteristics (e.g., domain wall energy density c,, 
and/or saturation magnetic moment M,). In this paper we use a varia- 
tional model to analyze the static stability of half-bubbles in materials 
having such gradients. Our model for the half-bubble and its instability 
modes is introduced in Section II. In Section III, we examine the stabil- 
ity of the model when the gradient in M/, is zero. While no stability 
occurs in this simple case, the calculation of Section III is a useful 
preliminary to the more general calculation of Section IV, where the 
gradient in M, is assumed to be non-zero. The results of the stability 
calculations and their interpretation are discussed in Section V. 





Fig. 1—Half-bubble magnetic domain configuration. 
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II. HALF-BUBBLE MODEL 


In constructing our model we idealize the half-bubble in two ways. 
First, since the details of the shape of the half-bubble are not known, 
we approximate the actual shape of the domain by a right circular 
cylinder of radius r, and height h (h assumed to be less than the platelet 
thickness), as shown in Fig. 2. This approximation gives the model 
domain a form identical to that of the standard cylindrical domain, the 
stability of which has been calculated by Thiele.* We shall exploit 
this identity in calculating the half-bubble energy and its derivatives. 
Second, we defer all consideration of the microscopic structure of the 
domain wall of the half-bubble and assume that the wall energy density 
o, varies linearly with distance through the platelet. 

The magnetostatic energy of our half-bubble model possesses an 
important invariance property. As we show in Appendix A for the 
case in which the M, gradient is zero, the magnetostatic energy of a 
cylindrical domain of fixed size located on the bottom surface of a 
platelet is independent of the position of the top surface of the platelet. 
This invariance is equally valid if M, is a function of z, or if the domain 
assumes a more general shape. Thus, if we conceptually move the upper 
platelet surface down to the top of the model half-bubble in Fig. 2, 
we see that the magnetostatic energy of the half-bubble is equal to 
that of a standard cylindrical domain of the same size. 

This invariance property has two important consequences. First, in 
the limit where the M, gradient is zero, it allows us to obtain the 
magnetostatic energy of our model half-bubble from the results of 
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Fig. 2—Model for the half-bubble magnetic domain. 
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Thiele® and Thiele et al.* for the standard cylindrical domain. Second, 
it indicates that the platelet thickness is not a fundamental length of 
the half-bubble as it is for the standard cylindrical domain. 

A convenient alternative standard of length, easily defined in terms 
of material characteristics, is provided by the material length at the 
bottom surface of the platelet (2 = 0): 


lL = o,(0)/47M? (0). (1) 


In defining the o, and M, gradient parameters y,, and yy, we then 
scale the variations in wall energy density and saturation magnetic 
moment to this material length: 


oy(Z) = v(0)(1 ae 12) ; (2a) 


M,@) = Myo(1 +“). (2b) 


Experimental observation’ indicates that there are three modes of 
instability which are of importance for the half-bubble: collapse, 
run-out, and “run-through’’. The general character of the collapse and 
run-out modes is already familiar from Thiele’s stability calculation® 
for the standard cylindrical domain; the half-bubble collapse is not 
purely radial, however, but involves axial collapse as well. The term 
“run-through” is used to describe the transition from a half-bubble to 
a standard cylindrical domain which contacts both surfaces of the 
platelet. 

As indicated in Fig. 2, the size of our model half-bubble is specified 
by the radius r, and the height h. In considering stability against run-out, 
we include an elliptic distortion r, and specify the radial boundary of 
the domain by 


r(~) = T. + T2 COS 2¢. (3) 


The model will be in equilibrium if the first derivatives of the total 
energy vanish. In the absence of in-plane anisotropy, 0H 7/drz = 0 at 
ro = 0 by symmetry (cf. Appendix B). The equilibrium conditions 
are then 

aE, aE, 


20) = 0. (4a,b) 





The stability of the resultant equilibrium is determined from the second 
derivatives, stability occurring when 
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2 2 2 2 
(SH). Gt)... ~ Gea). > 6) 
(Se). > 0, (6) 
and 
(SF), > ° © 


If the inequality in eq. (5) is replaced by the corresponding equality, 
then the half-bubble is on the verge of either collapse or run-through. 
If (0 7/dr3)eq. = 0, then the model is about to undergo run-out. 


III. HALF-BUBBLE STABILITY WITH CONSTANT VM 7 


We begin the actual calculation of the stability limits for the half- 
bubble by considering the limiting case of constant M,. This case 
lacks practical interest since we show that it leads to no region of 
stability for half-bubbles. In this limit, however, the problem loses much 
of the analytic complexity which characterizes the general case in which 
the gradient parameter yy, ~ 0. Our method of solution—treating the 
second derivatives of the energy as polynomials in h/l with coefficients 
depending only on 27r,/h-thus becomes clearly visible here. 

As noted in the preceding section, our half-bubble model is markedly 
similar in form to the model used by Thiele® and Thiele et al.* to treat 
the standard cylindrical domain. For this reason, in the present calcula- 
tion we adopt the notation of these papers wherever practical. 

The total energy of the half-bubble, measured relative to that of 
the uniformly-magnetized platelet, is 


Er = Ey, + AEy as AE x (8) 


where £,, is the wall energy, (including the energy of the top closure 
wall), AH, is the energy of interaction with the externally applied bias 
field, and AH, is the magnetostatic energy. Since the half-bubble is 
independent of the magnetic material above it, we can describe our 
model in terms of the unit step function wu as 


M = i,[1 — 2u(r(e) — rlu@uth — 2) (9) 


where 7;(¢) is as defined in eq. (3). 
The wall energy, consisting of integrals of c,,(z) over the side and top 
surfaces of the half-bubble, assumes the form 
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B= [ ae) de [ (rte + (4 a a +o.) | 22a. 0) 


From eq. (10) it is clear that different wall energy densities for side and 
top walls could easily be included in our model. While such a distinction 
has intuitive appeal, experimental evidence on this point is lacking. 
We thus omit this added degree of freedom in the present calculation. 

The energy of interaction with the external bias field H = 7,H, 
relative to that for the uniformly magnetized plate, is 





2r 
AE = -2 [ M-Hav = M.Hh | r() de. (11) 
h.—b. 0 


If the elliptical distortion r, is set equal to zero, H,, and AK, become 


E,, = 2ar sia (0)(1 aes Yeh) + or tou(0)(1 28 veh) (12a) 
= 2rr,ho,(h/2) + arzc,(h), (12b) 
AE, = 2nr-hM,H. (13) 


It should be noted that the interaction energy AH, is proportional to 
the domain volume mr?h. This characteristic is shared by the gradient 
term in the energy of the top wall [eq. (12)], leading us to view this 
term as arising from an effective bias field. The non-gradient part of 
the top wall energy lacks any h-dependence, being proportional to mr? . 
It may be thought of as coming from an h-dependent bias field analogous 
to that introduced by Liu’ in connection with sputtered thin films which 
are exchange-coupled to a magnetic bubble material. 

The relative magnetostatic energy is obtained from the general 
expression 


w=5 fav fe ay 2G (14) 


in the form 
ive) Qa i's) 
AE, = M? [ rar [ ae | r’ dr’ 
0 0 0 


[ dy'[kk’ — Ilo? — (e +h *] (18) 
where 


k(r, ro(e)) = 1 — Qu(n(e) — 1”), (16) 
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p =r +r” — 2rr' cos (vy — ¢’). (17) 
This integral is given in closed form by Thiele et al.* for the case r, = 0: 
AE y = — rh 24M I(x), (18) 
Iz) = _ 2x (0 + (1 + ike _ E( w ) _ x( a )} 
on x 1-2 l-2 
(19) 


where x = 2r,/h and K, E are complete elliptic integrals. 

It is now a straightforward matter to differentiate Hp with respect 
to r,, h, and r,. Using the definitions for the normalized bias field H 
and Thiele’s® wall-force function F 


AH = H/4rM, , (20) 


r= GP =F (+5) z)-1], 


we may conveniently write the equilibrium conditions dL;/dr, = 0, 
dE ,/dh = 0 in normalized form as 


sg ebb 
f(y, 4 teh) at (py eh) 4 ot — Pe) =0, 22a) 


(14 2h) 4 tome + dell — 810) + F@) = 0. 22) 
To show the biasing effect of the top wall we write eq. (22a) as 


1/2) fa an 1 1] _ F(x) = 0 (23) 


where /(h/2) and l(h) are the material lengths as measured at the middle 
and at the top of the domain, respectively, and are defined by the 
obvious generalization of eq. (1). Equation (23) should be compared 
with the force equation for the standard cylindrical domain [Ref. 2, 
eq. (69)] 


7 + «H — F(x) = 0. (24) 


Equations (22a) and (22b) may also be used to obtain expressions for 
the values of y,, and H required to obtain an equilibrium state charac- 
terized by specified values of I/h and zx: 


Yu = 25 (1 (1 = z) mi = 12) — 2F(zx) (25a) 
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ie 1 Bee x) _ o(2 + 4 )n@) ra (1 + 2)re), (25b) 


The stability or instability of a given equilibrium point is determined 
by examining the second derivatives of the energy. Calculating these 
derivatives and removing a factor of 8°’ = 4rh-4rM? we find 




















oe att tH SO (26a) 
Bicol Lalied) thera +32, 
B of = i. _ 5 I(a) + «F(x) — e ae ; (26c) 

earn gate) +g +3) 
+358 -41@)+27@)-542. 26a) 
Operations with these second derivatives are simplified if we define 
G(e) = — Ta) + aF@) - = EO. @7) 

Using the relation® 

= I(a) — F(x) = HS.(a) + 38,(0)) (28) 


we may write G(x) in terms of the stability functions S,(7) and S.() 
of Ref. 3: 
G(x) = 3x(S.(x) — S2(x)). (29) 
After consulting plots’ of S,(x) and S.(z), it is easily seen that G(x) = 0 
for x 2 0. 
If we use eqs. (25a, b) to eliminate y,,, H from eqs. (26a-d), we 
obtain expressions for the second derivatives of Hy at equilibrium in 
terms of //h and x only 


JE y 11 i) 2 
a(S are ) = Be (: = 3) T - Gm), (30a) 
GRP) es LA 1 
B (z a.) tag sche) ~ G2), (30b) 





ve.) pee on ee z) 1 
a( Oe ie = Bh i! 9 oe 5 G(x); (30c) 
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a ah (1 1) i 
al Ore Jeg. A \4 a a + a Gi). (80d) 
The derivatives are now in a form suitable for comparison with the 
stability conditions, eqs. (5-7). Inspection of eq. (30d) shows that 


B(0°E 7 /072)eq. > 0 for all positive values of l/h so that eq. (7) is satisfied 
for this system. Equation (6) is satisfied if 


Ll . 6 ae te 
p< OG +3) - oy 
Finally, we require 
| (Be) (222) _ (Bs y | 
B Or, Jeg. \ OW? Jeg. Or, Oh) eq. 
jae Nv " -£)-t( 1) 
-i()( oar ae g/e@ > 0 (32) 


which may be factored’ into l/h > 0 and 





2\ 1 
; > (2 + 3)ee)(1 —a- | (33) 
However, it is easy to show that eqs. (31), (83) and the condition l/h > 0 
cannot be simultaneously satisfied. Thus, if y,, = 0, the half-bubble 
model has no stability against collapse and/or run-through. 

While its result is negative, this stability calculation illustrates the 
method used in this paper for separating the //h- and d/h-dependence 
of the half-bubble model and for finding its stability region. We first 
use the equilibrium conditions to remove y, and H from the expressions 
for the second derivatives of H,. The resultant expressions are essen- 
tially polynomials in //h with coefficients depending only upon x = d/h 
[cf. the right-hand sides of eqs. (30a—d) and (82)]. For fixed x, the roots 
of these polynomials provide us with the limits of the stability region 
(if one exists). In the next section, we generalize this method to include 
the effects of non-zero ya, . 


IV. HALF-BUBBLE STABILITY WITH yy ~ 0 


We now generalize the equations of Section III to calculate the 
stability of the half-bubble model in the presence of a gradient in the 
saturation moment M,. The wall energy £,, is unchanged from the 
expression given in eq. (10) and the interaction energy AF, only gains 
a new factor to become 
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AEy = m1 + then f ri(o) do. (34) 
0 


The expressions for the magnetostatic energy, however, become sub- 
stantially more complex than those given in eqs. (15) through (19). 
The magnetization distribution of eq. (9) is now replaced by 


M = i.M (1 + mer — 2u(rs(y) — r)]u@)u(h — z). (85) 


The corresponding magnetostatic charge density is 
—V-M = —M,()[1 — 2u(r(y) — 7] 


x {300 — 6@ —h) + ue | — ae - h) += + culeuh = 2)|}. (36) 


The first two 6-functions in this expression are the familiar platelet 
surface charge terms which occur in the stability theory of the standard 
cylindrical domain [Ref. 3, eq. (19)]. The third 6-function gives the 
correction to the charge on the upper surface, while the product of 
step-functions describes the uniform volume charge resulting from the 
gradient Yar - 

The magnetostatic energy integral of eq. (14) has a quadratic de- 
pendence on the charge density —V-M. As a result, we have terms in 
AE ,, up to second order in yx, . For coon we write 

AE xy = EW + YM. vy - Vit i (37) 
where ES? is identical to the magnetostatic energy of the preceding 
section. The quantity y,:,Z$? has two parts. The first arises from the 
interaction of the ya, = O charge distribution with the uniform volume 
charge and is zero by symmetry. The second, representing the inter- 
action of the ys, = 0 charge distribution with the excess charge on the 
upper surface, is easily evaluated. We find 


wB? = Ul RY (38) 


The final term y3,HS? may be written as 
2 o 2a ie) 27 
uli, = 3 (vat) mx) | rar [ do | r’ dr’ i dg’ (kk’ — 1) 
0 0 0 0 
h 
x {o" — 2h [ dlp’ + (@ — hy} 
0 
h h : 
=e nf a | dz'[p” + (2 — nih. (39) 
0 0 


The p * term in the curly brackets corresponds to the self-energy of the 
excess charge on the upper platelet surface; the z-integral following it 
arises from the interaction of this excess charge with the uniform volume 
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charge, and the double-z-integral corresponds to the self-energy of the 
volume charge. After evaluating the z-integrals, we obtain the following 
integral expression for the second-order energy: 


2 ao 2a ee) 2a 
VE? = (x4) (0) par / i / r! dr" / do'(kk’ — 1) 
0 0 0 0 


- a fa h? 1 
x (2. = et HY, (40) 


Setting r. = 0, we find 
(2) 1 


7? = 4 (4) eem2ay (een 


{e — («+ 2°)F(x) — (3 — x) wh. (41) 


The derivatives of the total energy EL’, are now used to find the range 
of stability of the half-bubble model. Experimentally, this stable range 
would be mapped out by selecting materials characterized by given 
values of y,, and y, and then varying # to find the stable region in the 
variables d/l and h/l. The analytic form of our equations, however, 
precludes following such a physical approach in the present calculation. 
Instead, we use an extension of the method followed in Section III. 
We use the equilibrium conditions eqs. (4a, b) to eliminate y, and H 
from the second derivatives and then treat the second derivatives as 
polynomials in y = h/l, the coefficients of the polynomials being func- 
tions of x = d/h and yx. This manifestly unphysical approach gives 
the desired stable regions with a minimum of algebraic and numerical 
complexity. 

The equilibrium conditions corresponding to eqs. (4a, b) are con- 
veniently written in the following normalized forms: 


L(y + z) 4% Sd +2) + v(t + Yah) = (1 le tah (a) 





PS 1 (eet) |= — (1 + 2)F(a) + 2t(e) | = 0, (42a) 
; fe vel =e z) + = 1 Ala “F rah) + F@) - 3 7 L) 


+e we |e )- +1 | — load ~I(x) = 0. — (42b) 


These equations may be solved to find the appropriate values of y,, 
and # for an equilibrium state characterized by given values of yx , 2, 
and y = h/I. For brevity we write F(x) and I(x) as F and I respectively, 
and define a denominator D by 
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D= 7E + vu(1 as aa (43) 
We then find 


Vo = De — 2) + (yu — OF + 120 "Dy + yu(—10F + 2227'Dy? 
+22 —c+aort (Sta) |y 
+5 2 i= —-(1+2)F + (3 + x)r |} » (44a) 
po m{ ered) [ose dip 
+ aul (E+ 3) - B+ Srv 
- vil 3 (4 + (= —d+e +01) +304 ot |}. (44b) 


The denominator D in both of these expressions is always greater than 
zero and finite in the region of interest, and thus contributes no sin- 
gularities or zeros to the calculation. 

The second derivatives of EL’, are developed in much the same fashion 
as in Section III. They are cast in dimensionless form by extracting a 
factor of 8°’ = 4rh-4rM?(0), and y,, and 7 are then eliminated using 
eqs. (44a, b). The resultant expressions assume the form of polynomials 
in y after being multiplied by D. Because of the length of the expressions, 
we give only the final polynomial forms: 


Er) _ _(2,1\, (& _ 1) 
pe Se" ) = -(2+5) + (Se a /4 
~ va[ (4-3) F - (2 -5)r + rv 


< | a 
- da 


+C-beye-Gayy  w 
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baled), ~ 2049 - [For ¥e (+3) 

+ of (0-2) - G-s)e+ ES) 

tail +5(- 3) 
-G-F-s)r+G-2-3er 





ool Se), = 300-9) + [26-4 (- 349) 


dF 15 2 23 
+ va] —e(1 _ 2) aS + (8 ry — =p 3 re ry 


(2p (-5-2ey 
[E+ (e248 -0-ner Fh 
(45c) 
oa Se). G+) + [so + Bh 
rrf-(0-) E+ 8G - ar ey 
H(-f+ Br +(B+2- Bip 
+k (16 + 72) 2 1 -f47-2)@ 


: ae re Re ee : 
ate +i eri t (Ty 24h a alr|y 
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The problem of finding the ranges of parameters yielding stability 
of the half-bubble model now becomes one of numerical analysis. For 
given values of x and yy we find the zeros of the three polynomials 


(2B), - C0 
ve ( Ore JeqgN OR eq. \Oto Oh) og.) ’ 
2 2 
Dol Se) a bel Se) 

A stable region corresponds to a range of y beginning and ending with 
one of these zeros and having the property that all three of the poly- 
nomials are positive for values of y on the interior of the interval. 

The results of this procedure are shown for yx = —0.025, —0.05, 
and —0.1 in Figs. 3-5. The stable regions are plotted in the z, 1/y plane; 
1/y = l/h was used as a coordinate in order to facilitate comparison 
with the results of Thiele for the standard cylindrical domain (cf. Ref. 3, 
Fig. 3). In the interior of the stable region we have provided two families 
of curves which describe the behavior of a stable half-bubble. The 
families y, = constant and H = constant are obtained from the poly- 
nomials of eqs. (44a, b) after the boundaries of the stability region have 
been found. Since hf is not fixed in the half-bubble, the domains do not 
follow lines of 1/h = constant as A is varied. Rather, they follow the 
curves y,, = constant which characterize the given material. 











V. DISCUSSION 


We have explored the static stability of the half-bubble magnetic 
domain using a variational model. In spite of the simplicity of our 
model, the generality of some of the principles underlying our calcula- 
tion lead us to expect that our results have at least semi-quantitative 
validity (cf. Appendices A and B). 

Typical results of our stability calculation are shown in Figs. 3-5. 
The boundaries of each of the half-bubble stability regions shown in the 
figures roughly resemble two connected parabolic line segments, both 
parabolas opening to the right. The lower segment is determined by the 
zeros of the polynomial D8(0°E,/dr3).4, and corresponds to the run-out 
instability mode. The upper segment arises from the zeros of 


Es) (22) (22 ) | 

292 T Tv T 

me ( Or, Jeg. \ Oh” / ea. Or, Oh/ eq. 

This part of the boundary gives the limits of stability against run- 
through and collapse. Thus, the lines y,, = constant on the right side 
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Fig. 3—Static stability chart for the half-bubble model, yar = —0.025. 
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Fig. 4—Static stability chart for the half-bubble model, y;, = —0.050. 
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of the figure go from run-out at the low-H end of the line to collapse 
on the high- end. Values of y, which are algebraically small enough 
to correspond to lines on the left end of the stable region follow a different 
pattern, going from collapse to run-through as Af is decreased. Our 
stability charts then indicate two fundamentally different types of 
behavior as A is decreased: the half-bubbles may undergo run-through 
to assume the shape of a standard cylindrical domain (and then pre- 
sumably run out as A is decreased further),* or the half-bubble may run 
out first (presumably to undergo run-through as 7 is further lowered 
and give the usual demagnetized stripe pattern). These two types of 
behavior are separated by a critical wall energy gradient y,,, which is 
a function of ys, , being approximately —0.027 and —0.04 for yu = 
—0.025 and —0.050, respectively. (For y,;, = —0.1, the problem is 
more complicated and will be discussed below.) 

Comparison of the half-bubble stability regions with those of the 
standard cylindrical domain (Ref. 3, Fig. 3) reveals one marked dif- 
ference between the two systems. In the half-bubble case, there is a 
minimum stable value of d/h for given yy. For yx, = —0.025, —0.050, 
and —0.100, these minimum values are d/h = 0.39, 0.64, and 1.14, 
respectively. Long, narrow bubbles are then not attainable in the half- 
bubble system. 

It should be noted that all three of the stability charts correspond 
to yar < 0. While a general proof of the instability of half-bubbles with 
positive ys, has not been formulated, numerical investigations show no 
regions of stability for y;, > 0. As was shown to be the case for yx, = 0, 
eqs. (5) and (6) could not be satisfied simultaneously for yy, > 0. This 
property of gradient-stabilized half-bubbles is actually an advantage 
for device work as it implies a useful type of stability. A half-bubble 
originating on the lower surface of a platelet, after being allowed to 
undergo run-through, cannot be squeezed down to a half-bubble residing 
on the upper surface as 7 is increased again. The domain is thus forced 
to return to its original state, the half-bubble state on the opposite 
surface being unstable. 

In addition to the stability limits imposed by the derivatives of the 
energy function EH’, , we place a lower limit on l/h by requiring o,,(h) > 0 
and M,(h) > 0. While the motivation for the restriction ¢,,(h) > 0 is 
physically obvious, the reason for the requirement M,(h) > 0 stems 
more from the limitations of our model than from physical considerations. 
Magnetic materials may be compensated so that M,(z) goes linearly 
with z from positive to negative values. In the presence of a bias field, 
however, the magnetization in such a material would be aligned so that 
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its variation could no longer be described in terms of a linear variation 
with z. Such a compensated material would then not be susceptible to 
the present analysis and would presumably give rise to domains with 
different characteristics than those considered here. Taken together, 
the conditions o,(h) > 0 and M,(h) > 0 yield the restriction 


l/h > max (—Yyu , —Yar): (46) 


In practice, it may be necessary to require that //h be even greater than 
the lower limit implied by eq. (46) since o,,(z) > 0 must be satisfied not 
only at the top of the half-bubble (¢ = h), but throughout the region 
of the platelet above the half-bubble (2 > h). (By a simple extension 
of Appendix A, we need not require 1/,(z) > 0 in the region above the 
half-bubble.) These restrictions indicate that our calculation is not valid 
for part of the lower end of the stability region obtained for yy = 
—0.100 (cf. Fig. 5). As a consequence, the critical wall energy gradient 
we discussed above cannot be defined in this case. 

The normalized bias fields H required to stabilize half-bubbles are 
somewhat smaller than those needed for standard cylindrical domains 
as a result of the self-biasing effect provided by the upper closure wall 
(cf. Section III). Without going to unreasonably large values of y,, , 
cases are found which are stable with H = 0 or even H < 0 (H parallel 
to the magnetization inside the domain). In addition, bias margins for 
the half-bubble are of the order of those for the standard cylindrical 
domain, typically falling in the range 0.04 = AH = 0.08. Thus if 
4rM,,(0) is 200 gauss, bias field margins will be of the order 8 to 16 gauss. 

To demonstrate the use of the stability charts in obtaining numerical 
information, we consider the case y, = —0.100 and y,, = +0.025, for 
which a half-bubble is stable in zero bias field. The variations in a, 
and M, between the bottom and the top of the half-bubble are easily 
computed after obtaining a value of //h from the stability chart. In the 
present case, //h = 0.152 corresponds to zero bias. Evaluating y,h/l 
and yxh/l, we find that o,, increases by 16.4 percent and M, decreases 
by 65.8 percent between the bottom and the top of the half-bubble. 

The collapse and run-out fields for y,, = —0.100 and y,, = 0.025 are 
found by interpolating between the H = constant curves in Fig. 5. 
We find 7.01. & 0.017 and Ay... = —0.023 giving a bias margin of 
AH = 0.040. The size variations of the half-bubble over its stable range 
are best measured in terms of the z = 0 material length defined by 
eq. (1). For any operating point in the stable range, h/l = (l/h)~* can 
be found from the vertical coordinate of the point while d/l = (d/h)(I/h)~* 
can be found as the ratio of the horizontal and vertical coordinates. For 
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yu = —0.100 and y, = +0.025 we find (h/l).o1, = 5.38, (€/D eon. = 
10.85, (h/l) pun. = 7.27, and (d/l), = 80.91. Thus, while d increases 
by a factor of three between collapse and run-out, h increases by only 
35 percent over this same range. Since the curves H = constant and 
Y.» = constant are nearly parallel at the collapse end of the stable range, 
much of this change in h occurs with a very-small change in A near 
collapse. Thus, the relative variation in h may be decreased considerably 
by sacrificing a small fraction of the bias margin. The relative smallness 
of this variation in h is particularly important if half-bubbles on the 
lower platelet surface are to be propagated by permalloy circuits at the 
upper surface, as it is to be expected that there is an optimum circuit- 
domain separation for half-bubbles much as there is for standard 
cylindrical domains.” 

In conclusion, we have shown that the half-bubble has a number of 
properties which are attractive for device applications. Its stability in 
the absence of external biasing represents a considerable advantage over 
the externally biased standard cylindrical domain. While some com- 
plexity is added by the requirement of gradient materials, the materials 
problem is simplified on another front by the fact that material thickness 
is unimportant for half-bubbles. Indeed, the half-bubble stabilizes on 
only one surface; if propagated from the same surface, irregularities in 
the second surface are of no consequence. These attractive properties 
of the half-bubble and its possible utility in avoiding the undesirable 
properties of hard-bubbles may more than compensate for the added 
difficulty in growing half-bubble materials. 
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APPENDIX A 


Magnetostatic Energy Invariance 


Demagnetization fields equivalent to those of the standard cylindrical 
domain may be conceptually created through the use of two neutral 
subsystems of charges. The first consists of two parallel planes with 
magnetostatic charge density +//7, and the second consists of two 
disks with charge density -2M, placed on a common axis just inside 
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Fig. 6—Source-charge distribution for the demagnetization field of: (a) The stand- 
ard cylindrical domain. (b) The half-bubble model. 


the first subsystem (Fig. 6a). The distances between the upper plane 
and the upper disk, and between the lower plane and the lower disk 
are allowed to approach zero. The demagnetization field of the cylin- 
drical domain is the superposition of the fields of the two isolated 
subsystems: 


H, = Ha + Hus (47) 


The magnetostatic energy of the domain (including the energy of the 
uniformly magnetized platelet) is the sum of the self-energies of sub- 
systems 1 and 2 and their mutual interaction energy: 


ies | H2 dv (48a) 
Sr 


7 x [ Hav a ~ | He ave & | BaF dV. (48b) 
Sr Sr 4r 


The transition from the magnetostatic field configuration of the 
standard cylindrical domain to that of the half-bubble model is achieved 
by moving the upper plane of charge upward (Fig. 6b). The field strength 
of subsystem one is unaltered by this change, while the self-energy of the 
subsystem increases linearly with the distance moved. The self-energy 
of subsystem two is, of course, unchanged, but remarkably the mutual 
interaction energy of the two subsystems is unchanged also. The latter 
invariance follows from the invariance of the field strength of subsystem 
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one and of the charge-neutrality of subsystem two (so that all field 
lines originating in the isolated subsystem two also terminate there). 
From these two properties it is easily seen that any lines of Hz. extending 
above the upper disk give a zero contribution to the mutual interaction 
energy of the two subsystems. If we now observe that the magnetostatic 
energy of introducing a domain into a uniformly magnetized platelet 
is just the sum of the self-energy of subsystem two and the mutual 
interaction energy, we see that this differential magnetostatic energy is 
independent of the position of the upper surface of the platelet (provided 
only that the domain remains inside of the platelet). 

It is easily seen that this proof may be generalized to the cases in 
which M, is a function of z or in which the domain assumes a more 
general form. 


APPENDIX B 


Energy of Half-Bubble Under Elliptical Distortions 


The condition (dH ,/dr2),,-. = 0 in the absence of in-plane anisotropy 
is satisfied not only by our simple half-bubble model but by any convex 
domain having radial symmetry at equilibrium. For example, consider 
the generalized half-bubble shape 


ro(y, 2) = 7.(z) + ra(z) cos 2 (49) 
where r2(z) < r,(z) and r,(h) = 0 (i.e., the half-bubble is closed at the 
top). 

If there is no in-plane anisotropy, then 
Ey(r2) = Ey(—712) (50) 


since the shape function r,(y, z) = 7,(z) — 7r2(z) cos 2g transforms into 
that of eq. (49) under rotation. From eq. (50) it follows that 


Er) 7 
( or)» = (51) 


Consequently, cross-derivatives such as (0°H,/dr. dh),,-9 are also 
zero for the general domain shape given in eq. (49). 
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Earthquake Environment for Physical 
Design: A Statistical Analysis 


By S. C. LIU and L. W. FAGEL 
(Manuscript received April 14, 1972) 


A statistical analysts of nationwide seismic activities is made to determine 
the level and characteristics of earthquake motions to be expected in different 
geographical locations. Such information is needed to identify the severity 
of the earthquake threat to telephone facilities and to specify adequate 
design requirements. This paper describes the essential seismic environment 
including the expected earthquake-magnitude levels and the corresponding 
frequency of occurrence for different setsmic-risk zones, and the free-field 
and in-building motions that might be induced by earthquakes. Stochastic 
models are used to analyze the earthquake-occurrence statistics, and to 
describe the induced random ground motions. Historical earthquake data 
are used to verify the theory and to generate information useful for seismic 
design. 


I. INTRODUCTION 


Telephone communications are so vitally important to public health 
and safety that special efforts must be made to prevent disruption of 
these services by a destructive earthquake.’ For this reason, and in 
order to protect extremely expensive telephone plant investments, it 
is Important to incorporate earthquake-resistant design into telephone 
facilities. The urgency of safeguarding communications systems against 
earthquake hazards has become even more evident as a result of experi- 
ence gained. from the recent San Fernando earthquake.” At present, 
earthquake design loads for telephone plants are estimated through use 
of structural design procedures such as those outlined in the Uniform 
Building Code (UBC).* These procedures may be generally satisfactory 
for building design, but are not appropriate for determining seismic 
loads for equipment housed in central office buildings.* Consequently, 
earthquake-resistant guidelines should be established for the physical 
design of communications facilities, and these guidelines should be 


1957 


1958 THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1972 


based on realistic estimates of the seismic environments to which 
telephone plants may be exposed. These environments and the necessary 
theoretical background for defining them are presented in this paper. 
The dynamic behavior of telephone structures in these environments 
and design criteria are not considered. 

A statistical study approach is adopted here and historical earthquake 
data’ are used both to verify the theoretical analysis and to identify 
the earthquake environments of engineering significance. The objective 
in Section II is to derive two types of statistics: macrophenomena, which 
deal with the single highest value of random parameters that charac- 
terize the earthquake occurrence; and microphenomena, which describe 
the local earthquake motion. In other words, it is a study of the statis- 
tical distribution of the magnitude (or intensity) of the largest earth- 
quake in a given area, as well as the distribution of the largest amplitude 
of the nonstationary, earthquake-induced ground motion. In Section III, 
the theory is then applied to analyze available earthquake data and 
to generate information that is important from the viewpoint of seismic 
design. Geographic regions in this country that have had moderate- 
to-high levels of seismicity are selected (Fig. 1),° and the distribution 
of the intensity of the maximum annual earthquake for each region is 
then computed. Earthquake accelerograms along with the expected 
response spectra for given magnitude levels are presented which can 
be used for dynamic response analyses of structures. 


II. STATISTICAL ANALYSIS OF EARTHQUAKES 


Statistical techniques have been applied to many aspects of seismology 
and earthquake engineering; for example, in the data processing that 
involves estimating the travel time of seismic waves and locating their 
origin (epicenter), or in the analysis of seismic actions from a given area. 
In this section, emphasis is on this latter aspect, i.e., the statistical 
analysis of seismic action in a certain area from both the global and 
local standpoints. These two are distinguished by comparing the time 
scales they involve. The global (or macro-) analysis deals with the earth- 
quake occurrence over a long period of time measured in a multiple of 
years. In this case, the earthquake-occurrence times and magnitudes 
(energies) are treated as random sequences. The local (or ‘micro-) 
analysis deals with the free-field seismic wave motion itself, which 
generally lasts less than a minute. For each particular event, the local 
ground motion varies randomly with respect to time and constitutes 
a stochastic process. The objective of this section is to develop stochastic 


EARTHQUAKE ENVIRONMENT 1959 


Yves: 
Le 
GO wey O “iA ) 























ZONEO—[___—'|.:~ NODAMAGE 

ZONE 1—~/77777._—s MINOR DAMAGE; DISTANT EARTHQUAKES MAY CAUSE DAMAGE TO STRUC-— 
TURES WITH FUNDAMENTAL PERIODS GREATER THAN 1.0 SECOND; 
CORRESPONDS TO INTENSITES Z AND WI OF THE M. M. SCALE 

MODERATE DAMAGE; CORRESPONDS TO INTENSITY WII OF THE M. M. SCALE. 


4 MAJOR DAMAGE; CORRESPONDS TO INTENSITY MII AND HIGHER OF THE 
M. M. SCALE. 





THE PROBABLE FREQUENCY OF OCURRENCE OF DAMAGING EARTHQUAKES IN EACH ZONE WAS 
NOT CONSIDERED IN ASSIGNING RATINGS TO THE VARIOUS ZONES. 


Fig. 1—Earthquake risk map by the National Oceanic and Atmospheric Adminis- 
tration® and UBC.’ 


models that are suitable for seismic-risk study, use them to interpret 
existing earthquake data, and extract information that can be used to 
establish loading environments for telephone structures located any- 
where in the United States. 


2.1 Stochastic Model for Earthquake Occurrence 


Earthquakes can be considered to be a series of events randomly 
distributed on a real line (representing time), and the sequence of 
original times {t,} forms a point process.’ It is further assumed that 
the joint statistics of the respective number of shocks in any set of 
intervals are invariant under a translation of these intervals; this 
implies that {t,} is a stationary point process. The stationary point 
process generalizes certain aspects of renewal processes; in particular, 
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the interval lengths 7, = ¢ — t,-, between successive events need not 
be independently or identically distributed. 

The simplest stationary point process is the Poisson process. Intui- 
tively, the process {¢,} can be approximated as a Poisson process if it 
represents rare events.* More rigorously it requires that 7, be inde- 
pendently and identically distributed and follow a negative exponential 
function. Because of the simplicity and fairly general assumptions of 
the Poisson process, many studies have been made to test its fitness to 
the earthquake sequences.” ° Unfortunately, results of these tests in 
most cases indicate that the Poisson process is inadequate to explain 
the time distribution of low-magnitude shocks. The main deficiency of 
the simple Poisson model is its inability to describe the aftershocks 
which are often triggered by a large main shock. To account for the 
occurrence of aftershocks, the simple Poisson model for the main events 
is generalized to allow for the occurrence of more than one shock in a 
time unit, and a rate function g(t) defined for ¢ = O and normalized 
to unit area is introduced into the generalized model to describe the 
distribution of the total number of aftershocks.’ The function g(t) can 
be either a decaying exponential or an inverse-power relation depending 
on the magnitude, depth, and location of the main event.*'’® This model 
is known as the trigger or cluster model for earthquake processes. 

However, for most practical engineering purposes, the simple Poisson 
model for earthquakes appears to be adequate. In practice, an engineer 
is concerned with the earthquake risk of structures located in some 
specific geographic areas. The risk depends heavily on the statistics of 
large earthquakes of these areas, and the omission of small earthquakes 
or aftershock processes is relatively unimportant in terms of earthquake 
risk. For this reason, the simple Poisson process” combined with the 
extreme-value theory” is used in this study to derive the distribution 
function of the magnitude (or intensity, peak ground accleration, etc.) 
of the largest shock in a time interval. 

Let N(m), the expected number of earthquakes in an area per unit 
time whose magnitude M exceeds m, be given by Gutenberg and Richter’s 
familiar equation” 


logio N(m) = a — bm, m>0oO (1) 
or, equivalently, 


N(m) = ae™’™ (2) 


* Rare events may be justifiable if we consider only deep earthquakes or very 
large earthquakes.® 
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where a = exp (a ln 10), 8 = b In 10, and a and b are constants. Then 


P’(m) = prob (M > m) = No Sie" (3) 
and 
P(m) = prob (M S m) = 1 — e*", m> 0. (4) 


Therefore, the probability density function of the earthquake magnitude 
is given by 


p(m) = dP(m)/dm = Be*”, m>Q0O (5) 


which indicates that the magnitudes of earthquakes are independently 
and negative-exponentially distributed with a parameter 6. From 
eq. (5), the mean magnitude of all earthquakes of magnitude m > 0 
and its variance are, respectively, 1/8 and 1/8*. The mean return 
period, i.e., the mean interval in years between earthquakes having 
magnitude M > m, is 


T,, = WES = (") Jv. (6) 


From eq. (4) and under the assumption that the number of earth- 
earthquakes in an interval ¢ follows a Poisson distribution 


p(n, t) = (At)” exp (—d)/n! 


with an average rate ) > 0 (note that a = XZ), it can be shown that 
the magnitude of the largest shock, denoted by y, in the interval has 
the distribution 


FYy, t) 


prob (max m < y) = p(n, OPO) 


exp {—Ad[1 — (1 — e™)]} = exp (—Ate™), y>oO (7%) 


which is called the distribution of largest values of the first kind. From 
eq. (7), the probability density function of y in ¢ years becomes: 


f(y, t) = dF (y, t)/dy = Mp exp (—By — Me™), y>O (8) 


A concise yet detailed analysis of the above extreme value model and 
its application in earthquake statistics has been given by Epstein and 
Lomnitz.”” It is noted that the most significant parameters in the 
postulated earthquake occurrence model are \ and f. Based on historical 
earthquake magnitude data, the values of \ and @ can be estimated by 
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a least-squares regression procedure or by using Gumbel’s extreme-value 
plot, on which F(y, ¢) should asymptotically fit a straight line. This 
latter approach is used in Section III to analyze data from various 
active seismic regions in the country. 


2.2 Stochastic Model for Earthquake Ground Motion 


Attention is now directed to the microanalysis of earthquake statistics; 
a description of a stochastic model which can be used for this analysis 
follows. The objective is to derive the probability distribution function 
of the single highest amplitude of the earthquake ground motion having 
a given magnitude, as well as statistics of its induced structural response. 


2.2.1 Multiplicative Process Model 


Various stationary and nonstationary models have been proposed to 
describe the earthquake ground motion.” The approach is to adopt 
here the multiplicative process representative for the ground accelera- 
tion x(t) given by 


a(t) = o(i)n(d), (9) 


where ¢(t) is a deterministic envelope function and n(é) is a Gaussian 
stationary process with power spectral density S,(w) and autocorrelation 
function R,(r). The following expression taken from Jennings, et al.”* 
represents the envelope function (see Fig. 2): 


a ee a 
g(t) = 91 ih Stst, (10) 
| eats ty < t < i : 


where u > 0 is a constant and é; is the total duration of the accelerogram. 
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Fig. 2—Envelope function for earthquake accelerogram. 
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Under the assumption that the process n(é) is the response of a linear 
ground filter with characteristic constants w, (representing natural 
frequency) and é, (representing damping) to a Gaussian white noise 
with uniform power spectral density S, it can be shown” that the 
expected number of overshoots of x(é) per unit time over the double 
barrier + = +a is 


1 a 2;2 
Myai(E) = = TO,b° exp (-33 26 2b sfs)| at exp (-2¢) 


wT ad 
- 5 ad erfc (4 mr )]. (11) 


where o, and o, are standard deviations of n(é) and 2(t) = dn/dt, 
respectively. 

It is also possible to derive the distribution F(z, t) of z, the largest 
value of x(t) in a duration ¢t, by again assuming that the probability 
distribution for 7 crossings in ¢ is Poisson with a nonhomogeneous rate 
Niu) = fo ma(r) dr. By this assumption, which is good for high 
crossing levels, it can be easily found that F(z, t) = exp [—N,,,(4)] and 
the density function is f(z, t) = dF/dt. It should be noted that the 
analytical solution to the probability density function f(z, t) is extremely 
complex even when ¢(¢) has a very simple expression. However, the 
mean and mean-square values of z in an interval ¢ and its variance can 
be evaluated numerically without much difficulty. Note that these 
quantities are time-dependent and are functions of S, the spectral 
density of the white noise that produces n(). Some useful approximate 
solutions of the above statistical quantities are given in Appendix A. 





2.3 Distribution of the Maximum Structural Response 


The effect of the above earthquake process on structures can be 
represented by its response spectra, the maximum responses of single- 
degree-of-freedom systems with natural frequency w, and damping &, . 
The primary interest is to find the expected response spectra and the 
associated variances. Since the input process x(¢) is nonstationary, the 
response process u(t) is also nonstationary. The exact statistical behavior 
of u(t) is, therefore, very difficult to obtain. However, results of practical 
value can reasonably be obtained by considering either a stationary 
approximation or a Monte Carlo approach. The latter is more accurate, 
of course, and is used in this study; however, the stationary approxima- 
tion method is sometimes preferred because it can provide quicker and 
less-expensive numerical results. 
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2.3.1 Stationary Approximation 


This method provides the approximation of the expected response 
spectrum to an input earthquake process and follows directly from the 
same approximate treatment as that used in arriving at eqs. (13) to (17) 
of Appendix A. Here it is assumed that the initial build-up and final 
decaying portions of z(t) do not significantly affect the structural 
response u(é), and therefore, u(¢) can also be treated as stationary. 
On this basis, eqs. (13) to (17) will also be valid to approximate the 
corresponding quantities for the response process after replacing z by 
max |u|, o, by o., which is the standard deviation of u(t), and », by 
vy, = o;/(o,) = the zero crossing rate of u(t) from both above and 
below. The values of o, and o,; (the standard deviation of w(t) = du/dt) 
can be easily determined by integrating the power spectral density func- 
tions S,(w) and S,(w) = w’S,(w) respectively. The function S,(w) is 
given by S,(w) | H,(w) |’, where H,(w) = (w? — w” + 2it,w.w)* is the 
transfer function of the structure and based on the ground motion model, 
S,(w) = S/[(w” — w?)? + 4£,w7w?]. The results for the statistics of u(t) 
in terms of given input and system parameters are given in Appendix B. 


2.3.2 Monte Carlo Computation 


The second approach to establish the distribution of the maximum 
structure response is Monte Carlo computation; the procedure is 
demonstrated below. For each magnitude level M/, the expected spectral 
intensity SI). = fo? S,(Z. , 0.2) dT, , representing the damage (re- 
sponse) potential of an earthquake,” is obtained from Fig. 3, in which 
0.2 is the associated damping ratio, 7, = 27/w, = the natural period 
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Fig. 3—SI .2 Intensity and duration of strong shaking versus Richter magnitude. 
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of the structure, and the duration of the strong shocking ¢, is also given 
in this figure. The desired earthquake process x(é) is obtained by first 
generating the stationary process n(t) on a computer by the standard 
method,”* then shaping n(t) by $(é), and then normalizing it to match 
the expected spectral intensity SJ,... The response spectra S,(w, , &,) 
of each of the sample functions of x(¢) are then computed and averaged 
to obtain E[S,], where # denotes the ensemble average, and o,, , the 
corresponding standard deviation. By this procedure the numerical 
solutions converge to the true ones rapidly as the sample size increases. 
With efficient algorithms and computer facilities available, this method 
should prove satisfactory for practical purposes. This method is used in 
the next section to generate the desirable ensemble of earthquake 
accelerograms and the associated response spectra. 


III, APPLICATIONS 


The statistics of earthquakes that can occur within various specified 
geographic regions in the United States are next estimated using 
historical data and the theories given in the previous section. 


3.1 Return Periods of Karthquakes 


The distribution of the magnitude of the largest annual earthquakes 
that occur within a particular region is given by eq. (7) with ¢ = 1 year. 
Gumbel’s probability paper’ for extreme-value distributions of the first 
kind is used to plot this data. If the variate m does indeed follow a 
Poisson distribution, the distribution F(y, é) in eq. (7) will plot in 
Gumbel’s paper as a straight line. Earthquake return periods are equal 
to (1 — F(y, t))~’; therefore, the required information is automatically 
generated when the Gumbel probability paper is used. 

The procedure for constructing an extreme-value plot for a specified 
region is as follows. The highest annual modified Mercalli Intensity (2) 
[y in eq. (7)] for n years, as given in the U. S. earthquake catalogue,’ are 
tabulated in order of increasing size: y(1) y(2) S --- S y(n). For 
each y(z), the associated value of F(y(z), 1) = 7/(n + 1). The computed 
F(y(2), 1)’s are plotted on an extreme-value graph versus Richter 
Magnitude (/) by using the relationship M = 1 + 2I//3 as proposed 
by Gutenberg and Richter.” 

The geographic regions that have been referred to are selected after 
first dividing a map of the United States into 1-degree-longitude by 
1-degree-latitude rectangular areas (each segment is approximately 50 
by 70 miles) in Fig. 4. Each rectangle contains two numbers that are 
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Fig. 4—Earthquake intensity data map of the United States. 
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obtained from historical data’ if an earthquake of Modified Mercalli 
Intensity (J) greater than V had ever occurred within its boundaries. 
The arabic number represents the number of years when earthquakes 
have occurred, and the Roman numeral represents the largest intensity 
associated with any of those events. This information complements the 
seismic risk map in Fig. 1, and the two are considered together for 
defining boundaries around areas that appear to have comparable 
seismic history. These areas and boundaries are shown in Fig. 5a, and 
graphs showing Richter magnitude versus estimated return periods for 
earthquakes that may occur within these boundaries are shown in 
Fig. 5b. The largest magnitude distributions are plotted in Gumbel’s 
paper as straight lines in Fig. 5b with its slope being proportional to the 
standard deviation of the maximum magnitude levels. Note that these 
distributions are based on all earthquakes that occurred in the identified 
areas with the exception of the 1811-1812 New Madrid, Mo. and the 
1886 Charleston, South Carolina earthquakes. The return periods for 
these two violent events are estimated to be thousands of years; conse- 
quently, these are not regarded as meaningful data. Notice also that 
areas A, EK, and H (see Fig. 5a), although geographically discontinuous, 
have very similar and comparable seismicity (Fig. 4) and the calculated 
distribution data for these areas are represented by a common straight 
line in the Gumbel’s probability paper (Fig. 5b). Although various 
methods can be used to estimate the parameters of extreme distributions 
with different degrees of significance, such efforts are not warranted 
because a slight variation in the estimation is not an overriding concern 
in engineering risk analysis. Based on the same reason, distributions 
for areas B, F, G, I, J are also approximated by a common straight line. 

Based on the results of Fig. 5b, maps such as that in Fig. 6 showing 
the Richter magnitude levels associated with different return periods 
can be constructed. Such maps are of value because they not only show 
the different degree of severity of earthquake threat in various regions, 
but also give the information as to how frequently the damaging 
earthquakes would be expected to occur in these regions. It may be noted 
that similarities exist between Figs. 1 and 6 because the seismic-risk 
map was referenced in selecting the boundaries in the return-period 
study. Some of the dissimilarities are noteworthy, however: 


(c) The boundaries of the zone-3 area in parts of California and 
Nevada are not the same as those around the high-Richter- 
magnitude area in these states. The Great Central Valley below 
the 40th parallel has had greater seismic activity than most 
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Fig. 5a—Seismic regions and boundaries for magnitude distribution analysis. 
Fig. 5>—Earthquake Richter magnitude versus estimated return period. 


EARTHQUAKE ENVIRONMENT 1969 


y YY) 7 es 
: oy Yi) Vy X, 
V7 






















LY 


ec $e 






A y 





Fig. 6—Estimated maximum earthquake magnitudes for 40-year return period. 


other western states, zone-2 areas; therefore, its data is con- 
sidered with that of the neighboring zone-3 area. 

(iz) California and Nevada have earthquake histories that indicate 
that larger earthquakes will occur in this area than in other 
zone-3 areas for the same return periods. The estimated maxi- 
mum Richter magnitude is 8, whereas it is 6 for the Seattle area, 
6.7 for the Montana-Iowa-Utah zone-3 area, 6 for the St. Louis 
area, and 5.7 for the Boston and Charleston, 8. C. areas. Boston 
and Charleston do not appear to be more earthquake-prone 
than their respective surroundings which is in constrast to the 
information in Fig. 1. 

(z2t) Western New Mexico appears to have more earthquake potential 
than the rest of the Western United States, zone-2 region. 

(iv) A small area within the western part of Texas is designated as 
zone 2 solely because of the M = 6.4, 8/13/31 earthquake at 
Mt. Livermore. Insufficient data prevents an extreme value 
distribution to be established; therefore, the region is not 
designated to be different in Fig. 6 from its surroundings. 


Ground shaking that may occur in territories that are designated by 
an M = 6 design level in Fig. 6 (zone-1 in Fig. 1) are generally expected 
to be too small to be of engineering-design interest. Ground motion in 


1970 THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1972 


these areas could result from small-magnitude, local tremors or large- 
magnitude distant earthquakes; however, it does not appear to be a 
significant threat. The 17 = 5 design level shown for these areas in 
Fig. 6 should be regarded as an upper bound that is not necessarily 
applicable, particularly in geographic areas that have never experienced 
earthquakes (see Fig. 4). 


3.2 Free-Field Ground Motions 


It has not been feasible for building engineers or seismologists to 
consider geological conditions and then make precise, a priori predictions 
of ground motions that will occur at any particular point on the earth’s 
surface during an earthquake. This is so because most dependent 
parameters such as earthquake magnitude, focal depth, epicentral 
distance, slipped fault length, propagation path, and soil-media prop- 
erties, are usually not adequately known. 

A more direct and realistic approach to defining ground motions is to 
artificially synthesize them using characteristic parameters from past 
earthquake accelerogram measurements, such as dominant frequency, 
amplitude (Fig. 7), spectrum intensity, and duration of ground motion 
(Fig. 3), etc. Dominant frequencies of an earthquake are defined here 
as those corresponding to the dominant peak in the Fourier spectrum 
of the accelerogram. These are known to range between 6 and 60 
radians/s with a mean value of approximately 15 radians/s (Fig. 8). 
Peak acceleration amplitudes can be generally related to magnitude by 
the curve in Fig. 7. These are upper-bound values for compact alluvium. 
Accelerations on sites of different geology would fluctuate about the 
values shown in Fig. 7. 

The procedure and references for synthesizing an artificial earthquake 
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Fig. 7—Expected peak amplitudes of earthquake accelerograms. 
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Fig. 8—Typical fundamental frequencies of multistory buildings and dominant 
frequencies of strong-motion earthquakes. 


using the above parameters and a digital computer are given in Refs. 20 
and 23. A sample accelerogram for an earthquake of Richter magnitude 
6.3 is shown in Fig. 9a. The assumed values of the dependent parameters 
used to create this waveform are a dominant frequency of 15 radians/s, 
an SJ,,. of 2.3, and an envelope function which is comprised of a 3- 
second buildup, 9 seconds of strong stationary shaking, and an exponen- 
tial decay that lasts for 18 seconds. The waveform in Fig. 9a is one of 
an ensemble of 25 synthetic earthquakes generated according to the 
Monte Carlo method described earlier, and it corresponds to the average 
of the peak accelerations from the ensemble. This is 0.28 g’s which 
coincides with the acceleration value given by both Housner”™ and 
Gutenberg-Richter” for an earthquake of this magnitude. The velocity 
and displacement functions corresponding to this particular accelero- 
gram, obtained by integrations after performing baseline corrections, 
are shown in Figs. 9b and 9c. The peak amplitudes are 1.5 ft/s and 
1.3 ft respectively, 

Artificial earthquakes corresponding to different Richter magnitudes 
and dominant frequencies can be synthesized in the same manner. If the 
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Fig. 9—Acceleration, velocity, and displacement time-histories of synthetic 
earthquake with M = 6.3. 
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dominant frequency is close to the average (15 radians/s), it is simpler 
to approximate other horizontal ground motions by scaling the pre- 
viously mentioned waveforms with the factors shown in Table I which 
are the SI,.. values normalized to 1.0 for a magnitude 6.3 earthquake. 
This procedure will yield peak accelerations for different magnitudes of 
earthquakes that are approximately in agreement with Housner’s 
data” (Fig. 7). 

Figures 10a, 10b, and 10c are acceleration, velocity, and displacement 
response spectra of the accelerogram shown in Fig. 9a for damping 
ratios of 2 percent, 5 percent, and 10 percent. It is apparent that the 
acceleration, velocity, and displacement time functions of Figs. 9a, 
b and c together with these response spectra, are quantities associated 
with a particular sample member of the entire ensemble of earthquake 
ground motions. If a nondeterministic approach to structural analysis 
or design is used, the mean values and information about the variability 
of the ensemble response spectra are needed. The average response 
spectra and the corresponding standard deviations of the 25 synthesized 
earthquake samples are obtained by the Monte Carlo computation as 
described earlier and shown in Fig. 11. Note that these spectra represent 
only the elastic structural responses. For response analysis or design 
when inelastic behavior of structures is expected, the additional energy 
dissipation in the structural-foundation system should be taken into 
consideration. 

It is understood that the vertical motion of an earthquake is generally 
less severe than, and is not of as much concern to structural designers 
as the horizontal motion. The accelerograms and response spectra of 
vertical motions may be taken as one-half to two-thirds the horizontal 
accelerograms and spectra shown in Figs. 9 and 10, respectively. This 


TABLE I—Spectrrum INTENSITY SCALING Factors 


The following scaling factors are Spectrum Intensity Factors” (S7o.2) normalized 
to 1 for an earthquake of Richter magnitude 6.3. Free-field horizontal ground motions 
and response spectra may be estimated for design purposes by multiplying the ampli- 
tude coordinates on Figs. 9, 10, and 11 by these scaling factors. 


Richter Magnitude Scaling Factor 
0.6 


COMP crc 
NIo oN 

See OO 

Oe SOO 


1974 


DISPLACEMENT SPECTRA, Sg IN FT 


VELOCITY SPECTRA, Sy, IN FT/SEC 


ACCELERATION SPECTRA, Sq IN G’s 


THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1972 








0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 
NATURAL PERIOD IN SECONDS 


Fig. 10—Response spectra of earthquake accelerogram shown in Fig. 9. 
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Fig. 11—Expected velocity spectra and standard deviation of an earthquake 
process with M = 6.3. 


recommendation is made without the support of extensive theoretical 
analysis or data reduction as is done for the horizontal motion. However, 
it agrees with the limited data presently available, and it is consistent 
with the design practices adopted for other structures such as dams 
and nuclear power plants.” 


3.3 In-Building Motion Environments 


Telephone equipment installed within multistory buildings can 
generally be expected to encounter motion environments of greater 
intensity than equipment in single-story structures. Environments for 
the latter will be essentially free-field accelerograms, but multistory 
buildings amplify ground motions. It is not a trivial matter to specify 
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exactly how much amplification will take place, however, because the 
structural response of a building depends strongly on the properties 
of the structure and its founding soil. Soil properties affect not only the 
characteristics of the ground motions but also the fundamental frequency 
of the building. This in turn can drastically affect the way a building 
will vibrate and is generally referred to as soil-structure interaction.”°'*”” 

Typical fundamental frequencies for telephone buildings up to 
20 stories tall and founded on various soils are also given in Fig. 8. If 
one examines building-frequency formulas found in the UBC, it will 
be noted that some nontelephone buildings have lower frequencies 
than those shown in Fig. 8. Telephone buildings are usually extremely 
well designed and constructed, and in spite of the massive equipment 
contained within them, their natural frequencies are generally higher 
than those of conventional office buildings. 

It is important to understand that tall buildings, even the rigid ones 
constructed for telephone service, usually have fundamental frequencies 
that are below the dominant frequency contained in most earthquakes. 
This is fortunate because when the dominant frequency of an earthquake 
matches the fundamental frequency of a building, the acceleration 
responses of the structure can be much greater than in situations when 
these two frequencies do not match. This fact is illustrated in Fig. 12, 
which resulted from a study of the potential responses of a twenty-story 
building when subjected to earthquakes with various dominant fre- 
quencies.” It should be noted that the amplification of horizontal ground 
accelerations would clearly be highest when the earthquake’s dominant 
frequency happens to match the structure’s fundamental frequency and 
the building is founded on a stiff soil (V, = 4000 ft/s). This hypothetical 
situation is unlikely to occur however; in fact, historical data that are 
shown on the frequency coordinate of Fig. 8 suggest that a frequency 
match should not be expected for telephone buildings that are over 
eight stories. A frequency match for buildings under eight stories will 
not result in amplifications that are as large as those given in Fig. 12. 
The estimated maximum ratios of in-building acceleration to horizontal 
ground motion that can be used for design purposes is shown in Fig. 13, 
which is a composite plot applicable to all telephone buildings up to 
20 stories tall. These analytically-derived amplification factors may be 
regarded as upper-bound values that take into account the effects of 
the natural frequencies of the buildings, expected dominant frequencies 
of potential earthquakes, and all possible soil conditions on which 
buildings would conceivably be erected. Approximate in-building 
horizontal accelerograms can be artificially synthesized in a manner 
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Fig. 13—Estimated practical upper bound for in-building horizontal accelerations 
during an earthquake. 


similar to the generation of free-field motions, except that the value 
of SZI,.2 must be first multiplied by an amplification factor taken from 
Fig. 13. As an alternative to that procedure, the waveform shown in 
Fig. 9 and the corresponding response spectra in Figs. 10 and 11 may 
be multiplied by the appropriate scaling factor from Table I and also 
by the amplification factors given in Fig. 13. Note that these approxi- 
mations are valid assuming there is no resonance between the building 
and the floor-mounted equipment. 


IV. CONCLUDING REMARK 


An extensive statistical analysis of seismicity data and earthquake 
motions has been made. Simple yet realistic stochastic models are used 
to describe the earthquake occurrence process and the random local 
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ground motion. Information useful for seismic design of the telephone 
system is generated according to the theory presented and through 
the use of an available strong earthquake catalog. 

It should be pointed out that the main concern of this study is to 
provide information to assist in the physical design of earthquake- 
resistant structures by describing a realistic nationwide earthquake 
environment in meaningful engineering terms. The next step will be 
to establish appropriate structural design criteria. A simple approach 
for this would be to base such criteria on the most severe earthquake 
environment corresponding to the structure’s expected service life 
(say 40 to 50 years) and knowledge of its stress, strain and deformation 
tolerance limits under dynamic loads. For important structures such 
as telephone central offices which house various sensitive and expensive 
electronic communications equipment, a more rigorous approach based 
on both the seismic risk and cost analyses is desirable. An optimal 
design strategy in terms of earthquake and structure parameters can 
be reached by achieving a balance in the total construction cost and 
the expected loss due to earthquake damage. The mathematical formula- 
tion and detailed analysis of such an optimum seismic design procedure 
are reported in a forthcoming paper in the B.S.T.J.” 

Finally, it should be emphasized that the earthquake environments 
that are described are not intended, nor should they be construed, to 
be prophetic descriptions of future earthquakes. However, structures 
that are designed to adequately withstand these environments should 
consequently be expected to have a high probability of survival against 
earthquakes during their service life. 


APPENDIX A 


Approximate Solutions for the Peak Value Statvstics of x(t) 


Some approximate solutions can be obtained for the statistics of the 
peak amplitude of x(t) in the following manner. Assume that for ¢(é) 
in Fig. 2, the lengths of t, — 0 and ¢; — ¢, of the initial build-up and 
final decay are relatively short compared with the length ¢, — ¢, of the 
strong-motion phase of the accelerogram, and that the occurrence time 
t? of the extreme peak 2” of a sample function x’ (¢) will always occur 
in the range [t, , f.]|. These assumptions are justified from analyses of 
the behavior of structures under earthquake motion. Under these 
conditions and setting ¢ = 1, eq. (11) becomes 


Mya\(t) = v, EXP (-35) (12) 
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where », = o;/(to,) = w,/m is the expected number of zero crossings 
of the process n(é) per unit time from both above and below. The 
problem is then reduced to a stationary one. Letting tf = t, — t, and 
using the logic that follows eq. (11) readily lead to the following approx- 
imate solutions: 


Fe) = exp (—»,te*'/""") (13) 
n' t : —z? o 2 
ig) = a exp (-4 — pte? ) ; (14) 
zz(4+%)o,, (15) 
2’ = (A? + 2yo0r), (16) 
and 
wo. 


where A = (2 In » fo)’ and y is the Euler constant. If the ground 
acceleration is represented by a filtered white noise of constant power 
spectral density S, it can be shown that o, = (Sm/2E,w*)’” and that 
o;, = W,o,, Where w, and é, are the natural frequency and damping 
constant of the simple ground filter. From eqs. (15) and (17) and the 
expression for o, , 2 and o, can be expressed in terms of t, and S. 


APPENDIX B 


Stationary Approximation for the Statistics of Response Process u(t) 


The results for statistics of the response process u(¢) using the station- 
ary approximation method are: 


2 rSyB 


oo AaB aay ue 
22. wSyAs 
op = AA? a 
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al Avs) 
oS ae ( B (20) 


where A, = ww? , A, = 2w,w, (Ew, + &.), Ar = wy + wo? + 4E,E,w0, , 
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= 2(E.0, tp E Wg), B= A.A; -_ A, ? and 


s a 2E wee A 
Mo a(A? + 7)" 


The quantity Sy, in eq. (21) is the expected uniform power spectral 
density of the white noise which will produce an earthquake process 
a(t) of magnitude M with an expected peak ground acceleration Z,, 
[see eq. (15)]. 

Now the response spectrum can be expressed in terms of 2,,. The 
relationship between Z,;, and JM is shown in Fig. 3 based on the data 
from Housner.”® The expected (pseudo-) velocity spectrum, defined 
as S,(w, , &) = w, sup |u|, is then given by 

t 


(21) 


ES,(@, , £,) = wo, ? (22) 


where K = (2 In »,t,)' + y(2 In »,t,)7'”. The standard deviation of 
the pseudovelocity spectrum, in analogy to eq. (17), is given by 


g,, = 0 sup |u| = 1.28w,0,/K. (23) 
t 
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This paper ts concerned with the study of data multiplexing techniques 
which provide local distribution for user populations whose source char- 
acteristics may be categorized as Inquiry/Response. The techniques studied 
are Polling, Random Access, and a Loop System. In each method a group 
of users is multiplexed onto the same line which is connected to a Central 
Processor. The Central Processor forms the interface between the local 
system and long haul facilities, but does no computation on arriving messages 
other than to route them appropriately. The advantage of such systems is 
flexibility of operation and economy through sharing of equipment at 
the Central Processor. We focus our attention on the average round-trip 
message delay and use this as a measure for comparisons of the three 
techniques. 

‘Source models, in terms of calls per busy hour and number of bits per 
message for users as well as for computer responses, are developed which 
are appropriate in an Inquiry/Response context. Other factors taken 
into consideration are transmission rates, synchronization delay time, 
and allocation of available capacity to the transmission of overhead 7n- 
formation, such as polling messages and acknowledgment traffic. 

The results of the study are presented in graphic form, where the average 
delay due to traffic in the system is plotted as a function of the number of 
stations. Principal conclusions can be summarized as follows: 

Polling—It 1s found that this system is sensitive to the synchronization 
delay which takes place each time a user station transmits to the central 
facility. For reasonable choices of system parameters as many as 100 
stations can share the Central Processor without exceeding an average of 
1-second round trip message delay. 

Random Access—This system is not sensitive to synchronization delay. 
Over most of the range of load parameters, the Random Access system 
shows lower average delay than the other two systems. 

Loop System—The study revealed that this system has average delay 
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performance comparable to the other two. The Loop system is not sensitive 
to synchronization delay. 


I. INTRODUCTION 


One of the major functions of data communication networks is 
providing means by which users can access a distant computer in an 
interactive manner. Time-sharing, inquiry-response and credit checking 
are examples of applications which have a common communications 
requirement. They require rapid set up times for access to their re- 
spective computers and, as a consequence of their interactiveness, 
often require high-capacity return channels for return traffic. 

Providing access to a distant computer involves providing two 
major data transmission media, viz. long haul facilities and local distri- 
bution. (See Figs. 1 and 2.) While there are important problems as- 
sociated with each part, we shall focus our attention on local distri- 
bution techniques. We are motivated in this direction by the relatively 
high cost of providing local distribution compared to the total cost 
of transmission. 

The results presented in the sequel are most appropriate to the 
situation where users desiring access to distant computers are geo- 
graphically clustered. In this instance, traffic from several users is 
concentrated at what are designated as User Stations. We shall con- 
sider multiplexing techniques whereby traffic from User Stations is 
conveyed to a central point which we call the Central Processor. The 
role of the Central Processor is to route or direct message flow to and 
from User Stations and computers. 

This paper is devoted to an analysis of the roundtrip delay of three 
techniques which multiplex traffic onto common facilities. The analysis 
is based on models of user traffic and computer responses to user mes- 
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sages. In the applications we consider, traffic from the user is bursty, 
i.e., short periods of activity followed by relatively long idle periods. 
Each message from a user elicits a response from the computer. 

In the sequel we shall analyze the three alternatives for local distri- 
bution: 


(t) Polling—(See Fig. 1) Stations are polled by a Central Processor 
and transmit over a common line. Return traffic is multiplexed 
over a second common line by the Central Processor. 

(iz) Random Access—(See Fig. 1) Here User Stations transmit, 
at will on a common line, all messages that are generated. 
Positive acknowledgments are issued by the Central Processor 
upon error-free reception of the messages. Until an acknowl- 
edgment is received, the message is held in the station’s buffer. 
Again, return messages travel over a second common line. 

(iii) Loop System’—(See Fig. 2) User Stations share the same line. 
Traffic already on the line has priority, and newly generated 
messages are multiplexed into gaps in the line traffic. 


In each of these techniques, traffic to and from different User Stations 
share common facilities. The basic difference among the techniques 
lies in the means of controlling traffic accessing the system through 
the User Station. In the Polling System, the Central Processor exercises 
tight control over entering messages. As we shall see, in order to do this, 
overhead is incurred which causes an increased delay. In contrast, 
entering traffic in the Random Access System is loosely controlled 
by the Central Processor and delay attributable to overhead is reduced 
considerably. In the Loop System, as in the Random Access System, 
messages are multiplexed on the line at the User Station without direct 
control by the Central Processor. 

Our analysis of performance concentrates on buffering or queueing 
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delay. The results of this analysis will be used to calculate the round- 
trip queueing delay encountered by a user-generated message. A 
message going to the computer from a particular station is competing 
for the line with messages from other stations. Consequently it en- 
counters delay before it has sole access to shared facilities. Further, 
messages returning from the computer may be delayed in a queue 
at the Central Processor before being transmitted over the common 
line. We shall designate the sum of these two delays as the roundtrip 
delay. Of course, a message will encounter other delays such as propa- 
gation delay from the Central Processor to a distant computer and 
service time in the computer. However, these delays are independent 
of the local distribution system delays and can simply be added to 
the roundtrip delay that we calculate. 

A common thread running through each of these systems is inter- 
active queues. At each of the User Stations storage is assigned, if 
necessary, to queue up messages. However, since all the queues share 
the same server, the queues are not independent of one another. The 
exact treatment of interactive queues is mathematically difficult, 
therefore, in carrying out the analysis, certain approximations have 
been made. 

The three systems described above are compared on the basis of 
average roundtrip delay. Since we have computed the moment 
generating functions of forward and return delay, higher order moments 
can be found easily enough. From these higher order moments one 
can calculate other measures of performance. In so doing, the question 
of correlation between forward and roundtrip delay must be considered. 
We did not concern ourselves with this correlation since it has no 
bearing on the average roundtrip delay. 

In developing the models for the systems, we have attempted to 
take into account the constraints encountered in real communications 
systems. In the concluding section of the paper, examples using pa- 
rameters of the Digital Data System’ are presented. 


II. GENERAL CONCLUSIONS 


At the conclusion of the analytical investigations, we spend the 
final section of this paper in describing the delay performance of the 
systems. There also we compare, where possible, the advantages and 
disadvantages of various possible implementations. While the details 
appear later, it is worthwhile to give some general conclusions about 
the systems’ behavior here. 

For a given number of stations (V) and call rate (A), the average 
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message delay in the Polling System depends strongly on how quickly 
receivers can synchronize (s) as well as the available channel capacity. 
For instance, one can imagine that in the Polling System, polls take 
place over a separate channel of capacity 6. Now let us suppose that 
two separate systems could be designed, the first having zero synchro- 
nization time and the other taking 10 ms for the Central Processor 
to synchronize the polls. We would find as the number of stations 
increased, that the 10 ms spent at each station ultimately degrades 
the performance substantially. This effect can be seen by comparing 
Figs. 3 and 4 where we have depicted the above situation. For a large 
number of stations (~ 100) the two lower right families of curves 
rotate counterclockwise, indicating increasing delay. Note that even 
if we poll with capacity 6 = 100 bits/s (the upper right families of 
curves in Figs. 3 and 4), the delay is almost independent of the syn- 
chronization time and call rate. This is primarily due to the fact that 
we are forced to poll fewer stations (to maintain reasonable delay) 
and thus incur a smaller overhead induced delay. 

If the aim of the system is to serve a small number of stations each 
of which has a high calling rate, the synchronization time is relatively 
unimportant. Thus it may be desirable to use a low speed line with a 
small synchronization time. In Fig. 5, we depict this situation and 
note that it is possible in this case to achieve similar performance for 
a small number of stations (S$ 25) as is achieved with high speed line 
(Fig. 3). 

Since the Central Processor in the Random Access technique exercises 
loose control over traffic flow, there is a much smaller effect of over- 
head and correspondingly smaller delay (see Figs. 6 and 7). 

In the applications we consider, the volume of return traffic is greater 
than the volume of forward traffic, therefore, the first station in a 
Loop system suffers the largest delay. In considering the Loop system 
we focus on this station. Analytical results obtained in the sequel 
indicate that the Loop System compares favorably with each of the 
previous techniques. In the Loop System, such quantities as syn- 
chronization time and retransmission delay do not arise. What is 
relevant in this system is the processing delay, which involves examining 
the addresses of packets passing through the station. 


III. GLOSSARY OF IMPORTANT QUANTITIES 


N—number of User Stations 
}\—arrival rate of messages to the User Station (messages/s or 
calls/busy hour) 
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Fig. 3—Polling system delay in forward direction. 


By—message length generated at User Stations (bits) 
Cy—capacity of forward link (bits/s) 

m—time duration of message while on forward link (s) 
B,—return message length (bits) 

Cp—capacity of return link (bits/s) 
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M—time duration of message returning from the Central Processor 
(s) 
w—minimum time between polls in Polling System and mini- 
__ mum retransmission time in Random Access System (s) 
i, , /—mean and second moment of the time required to poll and 
readout messages in Polling System 
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Fig. 5—Polling system delay in forward direction. 


p—intensity of arriving traffic 
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dra—inbound message delay in Random Access System (s) 


d,—inbound message delay in Loop System (s) 
dp—inbound message delay in Polling System (s) 
drz—message delay on return link (s) 


é—capacity of return link allocated to polling messages in Polling 
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System and to positive acknowledgements in Random Access 
System (bits/s) 
s—synchronization time (s). 


Iv. SOURCE MODELS 


As indicated in the introduction we are primarily interested in the 
bursty, interactive user. Such users are encountered, for example, 
in credit checking and Inquiry-Response applications. The models 
that we shall develop are also applicable to time-sharing computer 
systems where users access a distant computer. Traffic characteristics 
of such users have been studied,** so that estimates of parameters 
characterizing the source (holding time, rate, etc.) are available. 

We model user traffic as consisting of Poisson arrivals of short, 
fixed-length messages. The arrival rate is designated as \ either in 
calls per busy hour or in messages per second. Each message consists 
of B, bits. If the line serving a User Station has a capacity of C, bits 
per second then m = B;/C,y seconds are required to transmit each 
message. 

In our analysis we shall assume that each User Station has unlimited 
buffer capacity so that buffer overflow never occurs. From the delay 
formulas that we develop, one can also derive buffer occupancy statistics, 
and thus imply the buffer size required for small overload probability. 
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Fig. 6—Random access system, delay in forward direction, fixed-time retransmission. 
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Fig. 7—Random access system, random retransmission, forward delay. 


For the message arrival rates that we consider, the station buffers 
need not be very large in terms of message lengths to achieve adequately 
small overflow probability. 

User messages arriving at the Central Processor are routed to a 
‘computer. Each message elicits a response from the computer. All of 
these return messages pass back again through the Central Processor. 
If messages arrive at the N user stations at average rate \ messages 
per second, then the average rate of return message flow through the 
Central Processor is NA messages per second. We assume that each 
return message has a fixed length of Bz bits, although our analysis 
can just as well be carried out for random message lengths. 

A basic assumption in our analysis is that the return message flow 
is Poisson, i.e., messages arrive at the Central Processor at random. 
Several conditions serve to randomize the return message flow. It 
may well be that messages sharing the same local facilities are destined 
for different computers, so that the returns from different computers 
will tend to be uncorrelated. The computers will be serving many 
local communities of users, introducing a random delay between suc- 
cessive messages returning to the same Central Processor. 

Under the assumption of a Poisson process of return messages, 
buffering must be provided at the Central Processor. Messages are 
read out of the buffer and transmitted to the User Station at a constant 
bit rate. 
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V. POLLING SYSTEM 


The polling technique is not new, and the reader is referred to 
Martin®’® for discussions on the present approach and others. The 
type of polling under study is ‘roll-call’ polling in which the Central 
Processor works down a list of users and, in sequence, interrogates 
User Stations sharing the same line. Upon being polled, the User 
Station may transmit messages. Since the generation of messages is 
random with respect to the time of polling, messages must queue in 
a buffer located at the User Station until a polling message arrives. 
As indicated earlier, we focus our attention on the average delay en- 
countered by a message. 

The arrangement of lines and User Stations for the Polling System 
is shown in Fig. 1. User Stations are bridged across common lines. 
Notice that, in contrast with the star connection, only one receiver 
and one transmitter is required at the Central Processor. The Polling 
System can be implemented using a “daisy chain” configuration in 
which the line between two user stations loops back to the Central 
Processor. 

Since polling messages are transmitted over a common line, User 
Stations must be capable of recognizing messages addressed to it. 
The same line is. used in common by messages returning from the 
computer. Again the User Station must be able to recognize return 
messages addressed to it. As we shall see, there is an interplay between 
forward and return message delay, in as much as polling messages 
and return traffic share the same physical line. 

The first component of the round-trip message delay that we con- 
sider is the delay at the User Station that a message encounters before 
it is put on the line. The calculation of this delay is complicated by 
the fact that the number of messages queueing at a particular station 
strongly depends upon the number of messages stored in all of the 
other stations. Since arrivals to any station are at a Poisson rate, 
the less frequently a station is polled, the more messages accumulate. 
Moreover, as messages accumulate, the polling cycle lengthens. 

The polling problem has been analyzed by Leibowitz’ by assuming 
a form of independence between queues.* The results of his analysis 
using this assumption compares well with the results of a rigorous 
analysis for a two-queue network. The nature of the independence 
assumption is such that the largest error occurs in a two-queue system. 

* A rigorous solution of the polling problem has recently been found by Eisen- 


berg.® However, in our application, which may involve a large number of stations, 
Leibowitz’s approach leads to more tractable results. 
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Of primary importance in the calculation of delay in the forward 
direction is a parameter that we designate as the ‘walk-time’ between 
stations. In the situation where a station has no message to transmit, 
the Central Processor must still spend time interrogating it. Time is 
required to transmit a polling message. The speed with which polling 
messages can be transmitted is limited, since messages returning from 
the computer share the same line. Further, after polling a station, the 
Central Processor must wait for a possible response. This time involves, 
for example, establishing synchronization between the User Station 
and the Central Processor. We shall return to the walk-time after 
we have considered its effect on delay. 

In calculating the delay of a message in this system, it is necessary 
to characterize the time between polls of a particular station. From 
the moment generating function derived in Appendix A we obtain 
the mean and the mean square values of the cycle time 


Nw 
-=T_Np (1a) 


_ NW = tw + pi)? + Nfw* + 2uiep + bmp), 
. 1—WNp 

In writing (la) and (1b), we take advantage of the fact that the walk- 
time is a constant value, w, and the message length is a constant value 
m = B,/Cy. Also we define p = \m. 

Messages that arrive at a User Station are buffered, awaiting the 
arrival of the poll. When the poll arrives, messages are read out of 
the buffer on a first come first served basis. Thus a message must wait 
until the station is polled as well as until messages that have arrived 
before it are read out of the buffer. In Appendix A, the generating 
function of delay is derived. From this generating function it can be 
shown that the average delay is given by 
t 
2. 
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dp = 





(1 + p) +m (2) 


where #, and # are given by eqs. (1a and b). 

Equations (1a) and (1b) show that the walk time, w, is an important 
parameter in the calculation of delay. We turn now to assigning values 
to it. Since the polling message must contain the address of the station 
being polled, it must be at least log, N bits long. Further, since it 
shares the line from the processor to the station with messages returning 
from the computer, not all of the return line capacity is available to 
transmit polling messages. We treat this situation by splitting the 
capacity of the return line. Suppose that of the Cz bit-per-second 
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capacity of the return line, 6 bits per second, interleaved with the 
message bits, are allocated entirely for polling messages. The minimum 
time required to transmit a polling message is then 


r = {[log, N] + 1}/6 (3) 


where [z] is the largest integer less than x. The remaining Cp — 6 
bit-per-second capacity of the return line is allocated to return message 
traffic. This allocation may be effected by either FDM or TDM. 

Other methods may be used to share the return line between polling 
and return messages. For example, messages may be transmitted over 
the same channel with priority given to polls or user traffic. The dif- 
ficulty with this approach is that there may be a large difference in 
the lengths of the two kinds of messages. Thus, even though a polling 
message is short and has priority, it may wait a relatively long time 
until a transmission of a return message is completed. 

Since User Stations continuously receive from the same point, the 
Central Processor, no time is required to establish synchronization 
for return or polling messages. This may not be true in the reverse 
direction since the Central Processor is receiving from a different 
station on each poll. Thus in the calculation of walk-time it may be 
necessary to allow time for the establishment of synchronization between 
the User Station and the Central Processor. Let this time be denoted s. 
The walk-time between stations is then the sum. 


w= s+ {[log, N] + 1}/6. (4) 


The effect achieved by splitting return line capacity on return message 
delay can now be calculated. Return messages are Bz bits in duration 
and arrive at the central processor at a Poisson rate of NX messages 
per second. The time required to transmit each message is Bg/(Cr — 5) 
seconds. The delay of a message in the Central Processor may be found 
from the analysis of an M/G/1 queue.’ It can be shown that the average 
return message delay is 


3, - —NABeMCe = OF 
¥ 2{1 — BprAN/(Ce — 8)] 


The average roundtrip delay is the sum of dp and dz . 


zi Br/(Cr 6). (5) 


VI. RANDOM ACCESS SYSTEM 


The Random Access system was suggested by the University of 
Hawaii’s ALOHA System.”° The configuration of this system is the 
same as shown for the Polling System in Fig. 1. As in the Polling System, 
messages arrive at each User Station at rate \ messages per second. 
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Messages are transmitted to the Central Processor on a first come 
first served basis. The length of a message in bits and the line rate are 
such that m seconds are required to transmit each constant length 
message. This constant length includes overhead bits required by the 
Random Access technique. 

Messages are transmitted at random by each station in the system. 
There is a probability of messages from different stations overlapping. 
In order to detect message overlap, parity check bits are transmitted 
along with each message and error detection is performed upon recep- 
tion. Only if a message is received error free is an acknowledgment 
transmitted to the User Station. Until a User Station recognizes a 
positive acknowledgment addressed to it, the message is held in a buffer. 
If, after a specified period of time, no acknowledgment is received, 
the message is retransmitted. This continues until an acknowledgment 
for a particular message is received. 

Messages returning from the computer are buffered and sent to the 
User Stations in the same manner as the Polling System. As in the 
Polling System, messages must contain addresses since messages 
destined for different stations share the same line. 

A basic relation governing traffic in the Random Access system has 
been derived by Abramson.’’’"* In this analysis, messages are assumed 
to enter the common line at Poisson rate. Suppose that a particular 
message is transmitted at time ¢t,. This message will encounter no 
interference if no other message is put on the line in the time interval 
t, + m. We have then 


Pr [message retransmission] = 1 — exp (—2Rm) (6) 


where F is the total transmission rate on the line, including retrans- 
missions. The rate of retransmission on the line is R(1 — e~7””). 

Now messages arrive at all user stations at rate \N messages per 
section. If there is no continuous buildup of messages at User Stations, 
we must have the relationship. 


R= N)\+ R(1 — exp (—2Rm)). 
This reduces to 
Nrdn = Rm exp (—2Rm). (7) 


An examination of eq. (7) discloses that the maximum information 
transmission rate is NAm = 1/2e & 0.18. As one attempts to exceed 
this, the retransmission rate increases so fast that less information 
gets through. 
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The basic assumption in this analysis is that line traffic, including 
retransmissions, forms a Poisson process. Simulation studies’” show 
that eq. (7) holds up to Nam & 0.1. For higher loadings, this simulation 
shows a much higher retransmission rate and, presumably, the Poisson 
assumption does not hold. 

A basie quantity in the calculation of delay is the time that the 
User Station waits before retransmission if no positive acknowledgment 
for a prior transmission is received. Let us consider how this ‘timeout’ 
period affects the system performance that a user might experience. 
If two messages do interfere, then each station must initiate a retrans- 
mission. If the times of each retransmission are the same, then inter- 
ference will persist. Thus some strategy must be used to avoid per- 
sistent interference as well as to provide for a short ‘timeout’ delay. 
To more fully appreciate what is meant by ‘short’ we have analyzed 
two strategies for retransmission. 

The first strategy assigns to each station a fixed or designated ‘time- 
out’ delay. This approach has the advantage that it completely avoids 
persistent interference. It has the disadvantage that some stations 
will have large (on the order of the number of stations times a fixed 
delay per station) delay. 

The second strategy (called randomized retransmission) attempts 
to make use of the fact that it is generally two stations that interfere, 
and that it is not necessary to distinguish between all NV stations in 
the system. This approach asks interfering stations to select a retrans- 
mission time by selecting it from a random sequence of retransmission 
times. If each station has a different sequence, then there is a small 
probability of persistent interference. The approach has the advantage 
of shortening the retransmission delay, but has the disadvantage that 
there is a non-zero probability of repeated interference. 

In the next two subsections we discuss the delay analysis for these 
strategies. Detailed calculations may be found in Appendix B. 


6.1 Fixed Timeout Delay 


Consider first the Random Access System implemented with fixed 
timeout retransmission. Suppose that a message is transmitted from 
a User Station at time ¢ = 0. Since we must allow for synchronization 
and for the transmission of an acknowledgment message, the message 
can be acknowledged no sooner than time t = s + {[log, N] + 1}/6 +m. 
If no acknowledgment is received at this time, station 7 waits for 
2im seconds (¢ = 1, 2, --- , N) and retransmits. This is repeated until 
a positive acknowledgment is received. If we assume that interferences 
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on successive retransmissions are independent, we have from eq. (6) 
Pr [message clearance time = w + m + kT; ] 
= Pr [k retransmissions] = [1 — exp (—2Rm)]‘ exp (—2Rm) (8) 
where w = {{[log. N] + 1}/6 + sand 7; = w + 2mzi. 


Average message delay at each station as well as higher moments 
of delay may be calculated by calling upon results from the analysis 
of an M/G/1 queue.” Messages arrive at a Poisson rate of \ messages 
per second and each message has a geometrically distributed service 
time as given by (8). From the generating function derived in Appendix 
B, it can be shown that the average delay for a message at station 7 is 


J _ 2[wP — 2Nm(1 — P)] — Nw*P + m* + QNm)'(1 — Py) 
os 2{1 — \NwP + m — 2Nm(L — P)]} 
(9) 
where P = exp (2Rm). 
The return message delay in the random access system is the same 
as in the Polling System. Fixed length messages of M bits return at 
rate N\ messages per second, with Cp — 6 bits per second capacity 


available for transmission. The resulting message delay is given by 
eq. (5). The roundtrip delay is the sum of dp, [eq. (9)] and dp [eq. (5)]. 


6.2 Randomized Retransmission Delay 


We turn now to consider a random retransmission technique. After 
an initial transmission, the User Station waits a minimum of s + 
([log. N] + 1)/6 seconds. If no acknowledgment has been received, 
the message is retransmitted after a randomly-distributed time interval. 
If no acknowledgment is received after this second transmission, the 
process is repeated. In our calculations we have taken the random 
timeout interval to be exponentially distributed with mean 1/a. 

The probability of interference on the initial transmission is given 
by eq. (6). On subsequent retransmissions the probability of inter- 
ference depends on a. For example, if a is very large, then retrans- 
mission for the two interfering stations occurs almost immediately 
after 


¢ eee 


and the probability of interference is high. As a decreases, the probability 
of interference decreases. We approximate the probability of inter- 
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ference on retransmission as 
Pr {retransmission interference] = 1 — exp[—2m(R + a@)]. (10) 
The average time between retransmissions is 
it, = s + {[log, N] + 1}/6 + l/a + m. (11) 


In Appendix B, the generating function of the time required to clear 
a message from a station’s buffer is derived. From this we can show 
that the mean and the mean square times to clear a message from a 
station buffer are respectively 





6b = w exp (—2Rm) + X(1/a + 2w) exp {+4m(R + a)} +m (12a) 
bv? = w’ exp (—2Rm) 
i" 2Xw* 2X(1 + 2wea)(1 + wa + waY) 
a exp {—4m(R + a)} a’ exp {—6m(R + a)} 
+ 2mb + m’ (12b) 
where 
X = [1 — exp (—2Rm)] exp {—2m(R + d)} 
and 


Y = 1 — exp {—2m(R + d)}. 


From the theory of the M/G/1 queue, it can be shown that the 
average delay is 
gab ee. 
o 2(1 — ab) 
where 6 and 8? are given by eqs. (12a) and (12b). 

Notice from (10) and (11) that decreasing a reduces the probability 
of interference while increasing the average retransmission delay. 
Thus there is an optimum value of @ that balances these effects, yielding 
a minimum average delay. 


+m (13) 


VII. LOOP SYSTEM 


The configuration of User Stations and Central Processor for the 
Loop System is shown in Fig. 2. As in the Polling and Random Access 
Systems, the line between User Stations can be looped back to the 
Central Processor. Traffic flow on the line is in terms of fixed size 
message slots. Messages arriving at a User Station are multiplexed 
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into these slots, one message to each slot. Traffic that is already on 
the line has priority; consequently a User Station must wait for an 
empty message slot. Traffic returning from the computer is addressed. 
At each station these addresses are examined. If a message is addressed 
to a particular station, it is taken off the line. Notice that in the Loop 
System, return and forward messages share the same line, consequently 
the total line length is half that required for Polling and Random Access 
Systems serving the same User Stations. 

Messages arrive at each User Station at a Poisson rate of \ messages 
per second. These messages are multiplexed on a first come first served 
basis. If a message arrives when the buffer is empty, it still must wait 
until the line is free before it can be multiplexed. 

Analytical and simulation results for this system have been presented 
in Refs. 18, 14, and 15. It was shown that the average message delay 
is approximated by the expression. 


7 Ey 0s) ee ee | 
t, = ml + pr) + 21 — p(l + pz)] = 1 211 + pr)[1 — p(l + px)] 


(14) 


where p, is the ratio of the average durations of line busy and idle 
periods. 7 and 7’ are respectively the mean and the mean square values 
of the durations of the line busy periods. In deriving eq. (14) we make 
use of the fact that the message length is a constant value of m seconds. 

A crucial quantity in the calculation of message delay is the line 
busy period. The characterization of the line busy period is simplified 
somewhat if we look at delay for what is usually the most critical 
station. In many applications the User Station receives more data 
than it transmits and as a consequence the traffic diminishes as one 
moves around the loop from the Central Processor. 

We have modeled the traffic to the Central Processor as consisting 
of Poisson arrivals of fixed length messages. Messages arrive at rate 
Nd messages per second and M = B/C, seconds are required to 
transmit each message. The busy period of the line out of the Central 
Processor is the busy period of an M/D/1 queue. It can be shown 
that” 





pa Pass (15a) 
1 — pr 
and 
2 
P= (B/C) (15b) 


7 dl = Pr) 
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where 
Pr = Brr»N/C2 . 
The average duration of a line idle period is 1/N\. We have then 


= BrNd/Cr _ Pr 


= ; 16: 
1 — pr 1 — pr re 


PL 
In order to provide a valid comparison with the other two systems, 
we must consider the processing time of a message at each station. 
A message, generated at the first station after the computer, passes 
through all of the other stations on the loop. At each station the address 
of the message is examined, entailing a delay of {[log. N] + 1}/C, 
where [2] is the largest integer less than x. Notice that since each station 
is receiving continuously from the same adjacent station, there is no 
synchronization delay. The cumulative delay of a message in going 
from the first User Station to the Central Processor is then 


d, =i, + N((loge N] + 1)/Cz (17) 


where #, is given by eq. (14). 

The return message delay for the Loop System is similar to the 
return message delay in the previous two systems. The delay in seconds 
is given by eq. (5) with 6 = 0. Thus the roundtrip delay is the sum of 
d;, [eq. (17)] and dp [eq. (5)]. 


VIII. EXAMPLES OF SYSTEM BEHAVIOR 


In this section, results of computations using the equations for 
average delay derived in the foregoing are presented. In presenting 
these results we were faced with the difficulty of choosing sets of pa- 
rameters that provide meaningful comparisons. There is such a wide 
latitude in the choice of parameters for each of the three systems that 
one could easily bury the reader in a mass of curves and tables. There- 
fore, we have limited ourselves to relatively few cases illustrating 
system behavior. It is not difficult to supplement the results we present 
here since the expressions we have derived for average delay are rela- 
tively easy to evaluate. 

Calculations have been made using values of user-related parameters 
(A, By and Br) which are appropriate in an Inquiry-Response context. 
Values for those parameters related directly to implementation (C; , 
Cz and s) are chosen with the Digital Data System’ in mind. 

The results of the computations are shown in the form of sets of 
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curves where average message delay is shown as a function of the 
number of user stations in the system. The rate at which messages 
arrive at User Stations, \, is given two values, 50 calls per busy hour 
and 500 calls per busy hour, representing, respectively, light and heavy 
calling rates. In the application we have in mind, messages in the 
forward direction tend to be short. We have taken this message length 
to be By = 250 bits. Three values were used for the lengths of messages 
returning from the computer 2.5 X 10°, 5 X 10°, and 10* bits. We note 
that this latter value is approximately the number of bits required 
to fill a CRT display. 

Since the systems we consider achieve economies by sharing, it is 
reasonable to put as many stations as possible in a system by choosing 
high line capacities Cy and Cz . In the Digital Data System, for example, 
56 kbits per second are available transmit to local stations over short 
distances. Thus we take Cp, = Cp = 56 X 10° bits per second. This 
choice is tempered by the fact that, in order to achieve this high rate, 
synchronous operation is required. As we have noted in connection 
with the Polling and Random Access systems, the Central Processor 
receives data from different stations on each transmission. A low 
estimate for the time required to adjust synchronization from reception 
to reception at this speed is s = 10 ms. As we shall see, this value of s 
may lead to high values of delay for the Polling System. Therefore, 
for comparison we examine Polling Systems where transmission in 
the forward direction is asynchronous with Cp, = 2400 bits per second 
and s = 1/2400 seconds. We have also examined fully synchronous 
operation where s = 0 and Cy = 56 X 10° bits per second. 

The parameter 6 comes into play in the Polling and the Random 
Access system. Recall that 6 is the portion of the return channel capacity 
allocated to the transmission of polling messages (Polling System) or 
positive acknowledgments (Random Access System). By varying this 
parameter, message delay in the forward direction is traded off against 
message delay in the reverse direction. For any particular system 
configuration, there is an optimum value for 6. However we have 
examined the effect of varying this parameter by choosing 6 = 10’, 
10° and 10* bits per second. 

Sets of curves of delay in the forward direction for the Polling System 
are shown in Figs. 3, 4 and 5. In Figs. 7 and 8, return delay is shown 
as a function of the number of stations. As one might expect, traffic 
characteristics are important in judging the merits of implementations 
of the Polling System. For } = 50 calls per busy hour, the 2,400 bps 
implementation (Fig. 5) yields lower delay than in the 56 kbps system 
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Fig. 8—Polling and random access system return delay. 


with 0.01-second synchronization time (Fig. 3) in almost all cases. 
In fact, at this calling rate the fully synchronous 56 kbps system (Fig. 4) 
has delay which is lower than that of the 2,400 bps system by a rela- 
tively small amount. This advantage disappears under the heavy 
calling rate (500 calls per busy hour). The lower speed implementation 
is far more sensitive to calling rates. 

Delay is sensitive to the parameter 6, which is that portion of the 
return channel allocated to polling messages. It is clear from that, 
for the traffie we consider, 6 = 10° bps is entirely too small. By in- 
creasing 6 to 10° bps, there is a large decrease in forward delay and a 
small increase in return delay. It is not clear that further advantage 
is obtained if 5 is increased to 10* bps. On Figs. 8 and 9, return delay 
may be very large for 6 = 10* bps and, for the same number of stations, 
low for 6 = 10° bps (e.g., N = 130 stations, MW = 2.5 X 10°/56 X 10° 
seconds on Fig. 9). For any given set of traffic characteristics there 
is an optimum value of 6 which minimizes total delay. 
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The results of the computation of average forward delay for the 
Random Access System are shown on Figs. 6 and 7. The curves of 
average return delay are the same as for the Polling System and are 
shown on Figs. 8 and 9. Unlike the Polling System, the effect of s and 
6 does not accumulate with the number of stations, and consequently 
is not very sensitive to these parameters. The Random Access System 
is more sensitive to calling rate than the Polling System primarily 
because calling rate affects the probability of a message being retrans- 
mitted [see eqs. (6) and (7)]. In Fig. 6, the results shown are for the 
station with the longest fixed timeout interval (¢ = JN) in eq. (11). 

The fixed-time retransmission implementation of the Random Access 
System (see Fig. 6) compares very well with the Polling System. For 
example, for s = 0.01 seconds, 6 = 10° bits per second and \ = 500 
calls per busy hour, the 100 station delay in the Polling System is 
nearly 0.9 seconds, (see Fig. 3) while the corresponding delay in the 
Random Access System is less than 0.2 seconds (see Fig. 6). The Random 
Access System with fixed-time retransmission also performs well in 
comparison with the s = 0 implementation of the Polling System. 

Recall that in the Random Access random retransmission strategy, 
the timeout interval is exponentially distributed with mean 1/a. By 
a process of trial and error we have found that a = Nd yields a rough 
minimum of average delay for VN > 10 stations. This value was used 
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Fig. 9—Polling and random access system return delay; \ = 500 calls/BH. 
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Fig. 10—Loop system forward and roundtrip delay. 


to obtain the results shown in Fig. 7. As we see from Figs. 6 and 7, 
this realization of the Random Access System yields performance that 
is superior to fixed-time retransmission. Random retransmission com- 
pares favorably with the best implementation of the Polling System 
(see Fig. 4). 

Our comparison may be biased somewhat in favor of the Random 
Access System since we have not taken into account overhead in that 
system. Recall that in order to detect errors, parity check bits along 
with information bits may be transmitted from the User Station. 
This will lengthen the message from the 250 bit message we have con- 
sidered. A larger value of m will cause increased delay by increasing 
the probability of retransmission [see eqs. (6) and (7)] and by increasing 
the timeout interval [see eq. (11)]. However, we felt that it would not 
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be necessary to lengthen the message very much, and consequently 
delay would be approximately the values we have shown. 

The results of average delay calculations for the Loop System are 
shown on Figs. 9 and 10. Forward and roundtrip delay for the station 
on the line immediately after the Central Processor are shown. Round- 
trip delay for the Loop System is even more sensitive to return message 
duration than the Polling or Random Access Systems. This is because 
messages going from the user terminal to the Central Processor are 
blocked by return messages. However, in the stable region below the 
knee of the delay curve, the Loop System gives performance com- 
parable to the best implementation of the Polling and Random Access 
Systems. For example, for = 2.5 X 10°/5.6 X 10* seconds, \ = 500 
calls per busy hour, and NV = 100 stations, the following average 
roundtrip delay estimates are obtained for each of the three systems: 
Loop System: 0.12 seconds (Fig. 11), Polling System: 0.16 seconds 
(Figs. 4 and 9), and Random Access System: 0.16 seconds (Figs. 7 
and 9). 
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Fig. 11—Loop system forward and roundtrip delay. 
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APPENDIX A 


In this appendix the moment generating function of forward delay 
in a Polling System is derived. Although other schedules may be ana- 
lyzed, we concentrate on the situation where each station is polled 
once and only once in a cycle. Messages arrive for multiplexing at the 
User Stations at a Poisson rate of \ messages per second. In carrying 
out the derivation, we shall denote the generating function of the 
message length in seconds as M(u). The generating function of the 
walk-time is denoted as W(u). In the text, we shall apply our results 
to the case where message length and walk times are constants. 

We consider first the ‘cycle time” of the polling sequence. This 
quantity is the time interval between polls of a particular station. 
As the polling sequence goes through a complete cycle, a random 
number of messages is encountered in each station’s buffer. The analysis 
is simplified considerably if we assume that this number of messages 
is independent and identically distributed from station to station. 
Under this assumption, Leibowitz’ shows that the moment generating 
function of the number of messages in each station’s buffer at the time 
of polling is given by 


P(x) = {W( — da) P@)}* (18) 


where @ = M(A — Az). 
Since messages arrive at a Poisson rate, a relationship between the 
number of messages in the buffer and the cycle time can be derived. 


Pr [n messages in a buffer at polling time cycle time = 7] 
= exp (—Ar)(Ar)"/nl. 
Averaging over the cycle time we have 


Pr [n messages in a buffer at polling time] 


_ i. exp (—dz)(dz)" 


nl p(t) dr 


where p,(r) is the probability density of the cycle time. By taking 
the Laplace-Stieltjes transform here it can be shown that, 
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P(x) = T,(\ — Ax) (19) 


where 7',(s) is the moment generating function of the cycle time. From 
(18) and (19) we have 


T.(u) = [(W(u)T.(s — M(u))]". (20) 


Differentiating (20) and setting wu = 0 yields the results on mean 
and mean square delay shown in eqs. (la) and (1b) respectively. 

Now a message arriving at a station must wait until the station is 
polled, and until all messages that have arrived before it have been 
multiplexed on the line. We first derive the generating function of the 
time the message must wait in the queue. We shall work under the 
assumption that messages arriving while prior messages are being 
multiplexed must wait until the next poll. Thus the queueing time 
of a customer can be written 


k 
dp = 7T + Se Mm; (21) 
i=l 
where 7 is the time interval until the next poll and m; ,7 = 1, 2,--- ,k 


are the lengths of k prior messages in the buffer. From (21) we can 
write the density of the delay dp as 


PY) = Pr[(T <dp S$ T+ dT] = i‘ dr Der — r)p(k, r) (22) 


where p,,(s) is the probability density of the message length and p(k, 7) 
is the joint density of k and 7. It can be shown from elementary prob- 
abilistic arguments that 
@o k 

p(k, r) = / aa, Pld Ne — a) exp {—Xa — 7)} (23) 
where p.(z) is the probability density of the cycle time. From (22) 
and (23) it can be shown that the moment generating function of the 
queueing delay is 


{A — AM(u)} — Tu) 


ie 
Trag(u) = i[u — \ + M(u)] (24) 


where 7,(u) is the moment generating function of the cycle time. 
In order to find the total message delay, we must add the multiplexing 
time to the queueing delay. Since these quantities are independent 
random variables, the moment generating function of the message 
delay is 
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Tr(u) = Trea(u)- Mu). (25) 


By differentiating 7’p(w) with respect to u and setting u = 0, the 
expression for delay given by eq. (2) of the text can be found. 


APPENDIX B 


In this appendix the generating function of message delay in the 
forward direction for the Random Access System is derived. Recall 
that after transmitting a message, a User Station waits for a positive 
acknowledgment from the Central Processor. The minimum time 
required to receive a positive acknowledgment is 


w = Hee AT ST ss, (26) 


If no acknowledgment is received, the original message is retrans- 
mitted. We consider two transmission strategies; a fixed transmission 
time which is different for each station, and a random transmission 
time. 

Under the fixed retransmission strategy, each station’s transmission 
interval is different by at most 2m seconds, thus the retransmission 
interval for the 7th station is 2mz. This interval obviates the possibility 
of the same two stations repeatedly interfering with one another. 
If & retransmissions of a message are necessary to clear a message 
from the station’s buffer, the total clearance time is 


r=wtm+ kw + mi) (27) 


for the 7th station. 
If message flow on the line is Poisson, then the probability of two 
messages interfering with one another is 


Pr [interference] = 1 — exp (—2Rm). (28) 


If we assume that subsequent retransmissions are independent trials, 
the total number of trials required to clear a message is geometrically 
distributed. We have 


Pr [k retransmissions] = [1 — exp (—2Rm)]* exp (—2Rm). (29) 


From (27) and (29) the moment generating function of the time re- 
quired to clear a message is 


S,:(u) = Elexp (—ulw + m + kT,])} = exp (—um)[exp (uw + 2Rm) 
— exp (2Rm — 2uim) + exp (—2uim)]"* (30) 
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where 
T; = WwW + 2m. 


In the language of queueing theory, S,,;(u) corresponds to the gen- 
erating function of the service time of a customer (message). Messages 
arrive for multiplexing at a Poisson rate of \ messages per second. 
The message delay can be found from the theory of the M/G/1 queue. 
We have 


(1 + AS7,(0))uS,(u) ‘ 


By substituting (30) into (31), differentiating with respect to wu and 
setting u = 0, the expression for the mean value of delay given in 
eq. (9) can be found. 

The derivation of the generating function of message delay for 
random retransmission proceeds along the similar lines. The retrans- 
mission interval is a geometrically-distributed random variable with 
mean value 1/a for each station. If k retransmissions are required, 
then the time to clear a message is 


Dralu) = (31) 


k 
r=(kK+Dwtmt+ Ye (32) 
j=1 
where é;,7 = 1, 2,--- , kare exponentially distributed random variables. 


The distribution of & is slightly different than in the case of fixed time- 
out retransmission. If two stations have interfered, there is a non-zero 
probability that they will interfere on the subsequent retransmission. 
We account for this phenomenon by approximating the probability 
of retransmission by the expression 


Pr [interference on retransmission] = 1 — exp (—2(R + a)m). (88) 
The probability of no retransmissions is 
Pr [no retransmission] = exp (—2Rm). (84) 
The probability of k retransmissions is 
Pr [k retransmissions] = [1 — exp (—2Rm)][1 — exp (—2(R + a)m)]*" 
‘exp (—2m(R + a)); &=1,2,---. (35) 


From (31), (33) and (34) it can be shown that the generating function 
of message clearance is 


DATA MULTIPLEXING TECHNIQUES 2011 


ll 


S2,(u) 


B| exp (—« + 1)w+m-— 5 a) 


exp (—un)| exp (—2Rm — wu) 


[1 — exp (—2Rm)] exp (—2m(R + a@)) exp (—200) |. 


a \ + u— dA exp (—wu)[1 — exp (—2m(R + a))] 


(36) 


The generating function of message delay can be found from the theory 
of the M/G/1 queue. The generating function of message delay is 
the same as in (31) with S,;(u) replaced by S2;(u). The first two moments 
of the time needed to clear a message given by (12a) and (12b) are 
found from successive differentiations of (86). The formula for average 
delay given in (13) is well known.° 
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On Fast Start-Up Data Communication 
Systems Using Pseudo-Random 
Training Sequences 


By R. W. CHANG and E. Y. HO 
(Manuscript received June 21, 1972) 


This paper analyzes the start-up performance of automatic transversal 
equalizers when mazximum-length pseudo-random sequences of short 
periods are selected as the training signals for fast start-up purposes. 
Single-sideband Nyquist systems are considered because they represent 
the limiting case of vestigial-sideband systems with small excess band- 
width. It 1s shown that the equalizer 1s capable of fast start-up except 
in some rare situations which can be avoided by using proper timing, 
phase, and equalizer initial settings. The results also show that the equal- 
ier tap convergence rate is independent of the phase characteristic of the 
communication channel and of the choice of the pseudo-random sequences 
which have the same period. 

The equalizer is set up in the training period by minimizing the mean- 
square error between the equalizer output and the transmitted pseudo- 
random sequence, which is different from the mean-square error for random 
data. Surprisingly, we have found that, even for pseudo-random sequences 
of very short periods, this start-up algorithm results in only a slight deg- 
radation in system performance. Accordingly, good system performance 
can be expected immediately after the system switches from the training 
mode to the data mode. 


I. INTRODUCTION 


Pseudo-random sequences have been used in the past as training 
signals for setting up automatic transversal equalizers during start-up 
periods.’’” For fast start-up, it is desirable to know how the equalizer 
settling time depends on the choice of the pseudo-random sequence; 
the channel characteristics, and the initial receiver conditions. These 
problems are examined in the first part of this paper for single-sideband 
Nyquist systems. We present basic theories from which the reader 
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can work out numerical examples of special interest. The results are 
compared with those obtained previously’ for a different class of training 
signals (isolated test pulses). An important difference between these 
two cases is pointed out. 

When pseudo-random sequences are used, it is most convenient 
to adjust the equalizer tap gains to minimize the mean-square error 
between the equalizer output and the transmitted pseudo-random 
sequence. It is not immediately clear how closely this simple algorithm 
optimizes the data set performance for transmission of random data 
(because an equalizer setting optimum for pseudo-random sequence 
transmission is not necessarily optimum for random data transmission, 
particularly when pseudo-random sequences with very short periods 
are used for fast start-up purpose). This problem is examined in Section 
IV and the analysis is illustrated by examples. 

Section V summarizes the results of this paper. The reader mainly 
interested in the conclusions and their implications may read Section 
V next. 


II. MATHEMATICAL MODEL AND FUNDAMENTALS 


An amplitude modulation data communication system with a con- 
ventional tapped delay line transversal equalizer is depicted in Fig. 1. 
During data transmission, the transmitter transmits the information 
digits, {d;}, sequentially at time instants ¢ = ---,4 — T, 4, 
t, + T, ---. The equalizer output is sampled sequentially at the symbol 
rate to recover the information digits. Let the 7th equalizer output 
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Fig. 1—Block diagram of an amplitude-modulation data communication system. 
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sample be y; . We adopt the familiar mean-square error (MSE) criterion 
and adjust the gain controls of the equalizer to minimize the MSE 
between y; and d; . We assume {d;} is an ergodic process, hence the 
mean-square error can be written as 


lim * 2, (y; — dy 
yoo i=1 
((y; — d,)”) (1) 
where (x) denotes the time average of x. 

It can be seen from Fig. 1 that 


& 


ll 


ui) = Yo Catt — (& — YN, (2a) 
and 
Hi = > d:h(t — 1), (2b) 


where A(t) is the overall system (without equalizer) impulse response. 
For the sake of simplicity, we shall shift the time origin and use the 
abbreviations y; = y(iT), 7; = x«(iT), and h; = A(iT). Thus (1) can 
be written as 


§ = C/AC — 2C'V + (d’), (3) 
where 
C, 
ct eel (4) 
Cy 
Qi1 Gio *°' Ain 
a ee a a (5) 
GQyi1 G@ne2 °*** Onn 
ai; Bopeitagy) 


lI 


2 > (din d,,)hes1—mh—ja1~n 


t= 1,2; tN, pH 12,028 N, (6) 
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Vi 

Ve 

Vy 
and 


Vie = (i-n41 di) 
be (di An)hi—eai—m ) k= Lop e984 N. (8) 


Let 0&/0C; be the partial derivative of & with respect to C; ,7 = 1 
to NV, and let 0&/dC represent an N X 1 column vector the 7th element 
of which is 0&/dC; . From (3), we obtain 

0& 

AC = 2AC — 2V. (9) 
The optimum value of C that minimizes the MSE & is denoted C,,, . 
It is clear from (9) that 


Cy Say: (10) 
Thus, the minimum value of 6 when C = C,,, is 
Emin = (03) — V’A'V. (11) 
Let e denote the difference between C and C,,, ; Le., 
=C-—A'V. (12) 


Now consider the adjustment of the equalizer. As is well known, 
the equalizer can be adjusted in the training period by transmitting 
either a succession of isolated test pulses or a special class of pseudo- 
random sequence (pseudo-noise sequence’). The case of sending isolated 
test pulses during the training period has been considered by Chang.* 
In this paper, we examine the case of sending pseudo-noise sequences." 

In the training period, a binary pseudo-noise sequence is applied 
to the transmitter input. Since an adjustment is made at the end of 
each training sequence period and since we wish to make the largest 
number of adjustments during a fixed training time, we consider the 
shortest* possible pseudo-random sequence period (i.e., the case where 

t The best known sequences of this type are the m-sequences (also known as 
maximal-length linear recurring sequences or maximum-length pseudo-random 
sequences). 


+ When the period of pseudo-random sequence is shorter than the length of the 
transversed equalizer, the analysis is difficult because A~! may not exist. 


EQUALIZER START-UP PERFORMANCE 2017 


the period of the pseudo-noise sequence is equal to the length of the 
transversal equalizer). 

Let the pseudo-noise sequence be denoted by 66283 --- By BiB283 -- 
GB. , where 8, is the last bit. From (2b), the input to the equalizer 
can be written as 


a(t) = » B,A(t — KT). (2b) 


Practically speaking, we may assume that A(é) is time limited. Then 
it can be shown that for ¢ larger than a certain value, say t) , x(t) will 
be a periodic function of period NT’; i.e., 


at) = a2t+ NT), t <t and t+ ANT Sh’T. (13) 


In the training period, the values of a;; , Vi , &, C,,, , A, and V are 
denoted by af , V* , 8&*, CX, , A*, and V*, respectively. From (6), 
we obtain 


* 
as; 


ll 


> > (so) | /Raratenees ear eran 
a ae hetat—-mh-j41-m+kn 


k 


N-1 


! 7 
~ N pe 2 ay heisi—mh—jai-meawer all 4, 7. (6) 


t=1 


From (8) 
Vi » (BmBi)Rimke1—m 


I 


1 N-1 
> hiusisin = N >; yy hiusisinet ’ k= i; 2, reg ,N 
7 l=1 


j = integers. (8) 


The partial derivative 08*/aC; , 7 = 1 to N, can be computed from 
each block of N samples of z(t), and the gain control C; is adjusted 
by an amount proportional to 08*/dC; at the end of each block. For 
example, 0&*/dC;; is computed from 7;,,; to 2:7 and C; is adjusted 
after r;.. . Then 08*/dC; is computed from 2,441 tO 2j4ey and C; is 
adjusted after x;..2» . The optimum tap setting that minimizes the 
MSE, &"*, is 


Cs = (A. v= 


We now proceed to examine the convergence of &*. 
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Ill. SINGLE-SIDEBAND NYQUIST SYSTEMS 


In this section, we consider single-sideband data communication 
systems which transmit at the Nyquist rate with sin x/z pulses (here- 
after referred to as single-sideband Nyquist system). Such systems 
are considered because they represent the limiting case of sharp rolloff 
vestigial-sideband systems. 

The transfer functions of the transmitting filter, transmission medium, 
and receiving filter (see Fig. 1) are F,(f), F.(f) and F3(f), respectively. 
The F’,(f) are of the following form 


F(f) = | F(f) | Cee 1=1, 2,3 (14) 


where J is used to denote the imaginary number V —1. 
In a single-sideband Nyquist system, | F',(f) | and | F'3(f) | are specified 
by 


I 
_ 


| Fi) | hRelflsh 
= 0 otherwise, (15) 
and 


I-FOl=1 ASI Sh 
0 otherwise. (16) 


In general, with lower single-sideband transmission, the carrier fre- 
quency, f, , is set equal to f, . Let H(f) denote the Fourier transform 
of A(t), which is the overall system impulse response at the equalizer 
input. It can be shown that 


A(f) = 3 | Pilf — fo) | 


eT IBL Fo) 48a Fo) +Ba S—Fe) “24 FF) FO Osfsh—hfi 

= 3 | F.(f + f.) | 

QTL 0) +Ba Fe) +B S4Fe) OR S4Fe) —(f, -f) sf <0 
2% otherwise, (17) 


where @ represents demodulating carrier phase. The signaling interval is 
: 1 

D ee ee 18 

2G — fi “ 

Since the time samples A(:T), 7 = --- 0, 1, 2, --- are taken at the 

Nyquist rate, we obtain from the sampling theorem, Parseval’s theorem, 
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and (17) 


gt oe )) 4 Rgepaaltayeiok 


m= 


2 — A) [ae — ar + TAG — + T) at 


et) [ (eos 2nfi — ATI FG - 1) af. 9) 


Substituting (19) into (6), we obtain 
N-1 

at = Dai —F+kN)-— Vg -F+EN +) all i,j. (20) 
k k t=1 


It is clear from (19) and (20) that a, is independent of the demodulating 
carrier phase 6, the system timing ¢, , and the phase characteristics 
B;(f) of the system. We also note that for a fixed NV, a is independent 
of the choice of the pseudo-noise sequence. [The pseudo-noise sequence 
8, , 82, °+°* does not appear in (19) or (20).] Using the method in Ref. 4, 
it is concluded that the equalizer tap convergence rate is independent 
of the demodulating carrier phase 0, the system timing ¢, , the phase 
characteristics of the system, and the choice of the pseudo-noise se- 
quence for fixed N. 

Note from (19) and (20) that a depends on the amplitude char- 
acteristic | F.(f) | of the transmission medium. Since amplitude dis- 
tortion is not severe in private line systems, in the following discussion 
we assume that 


IFAl=1 Asrsk. (21) 


Substituting (21) into (19), and neglecting a normalizing constant 
(fe — ae we obtain 


al a9 


= Q, a j. (22) 
Substituting (22) into (20) gives 
a¥, = 1, t=7 
1 , ; 
=-y ti. (23) 


Therefore, the eigenvalues of A* are 
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Malte, as ee ee | 
1 


The eigenvector uy corresponding to Ay is an N X 1 vector whose 
elements are all unity. Since all but the last eigenvalue are equal, 
the equalizer can settle rapidly except in the case where the initial &* 
contains a large component (ejuy)’Ay , where e; is the initial tap setting 
error vector®. Since \y is small and since (e/uy)” cannot be exceptionally 
large with proper timing, phase, and initial equalizer settings’, it is 
very unlikely that 8* would contain a large component (ejuy)*Aw . 
Therefore, the equalizer can settle rapidly in the training period. 


IV. FURTHER ANALYSIS OF SYSTEM PERFORMANCE 


At the end of the training period, the equalizer taps are set very 
nearly to C*, . The data set is then switched to the data transmission 
mode. Since statistics of the true data differ from those of the training 
pseudo-random sequence, the optimum tap settings, C*, , obtained 
for a training sequence cannot also be the optimum one for the true 
data. Thus, system performance degradation during the early stage 
of data transmission is expected, even if the data set is equipped with 
an adaptive equalizer. We now proceed to determine this degradation. 

We assume zero-mean independent information digits (binary or 
multilevel). The signal level is normalized such that 


di) = 1. (25) 
The mean-square error, &, can be obtained from (38), 
& = (Co.)ACH) — 2(C)/V + 1, (26) 


where A and V are given by (5) and (7), respectively. From (6), (8), 
(19), and (25), we obtain 


ai; = gv — 9), (27) 
and 
Vi, = henat . (28) 
The a* and V*% can be rewritten as 
az; = a; + A:,(N), (29) 


t For example, one may use the method described in Ref. 1. 
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and 
Vi = Vit 7.(N), (30) 


where 


N-1 


aN) = DoG- f+) — HD LoG-4i+mN+D, GY 


t=1 


and 


vlM) = Dhow — HD Dhan 
From (26), we have 
& = CiAC.. — 2C..V + 1 + (8C)’A(C) 
= &nin + (8C)’A(8C), (33) 


where C,,, is the optimum tap setting for the true data and 8C is the 
difference between C*, and C,,, , 


aC = Cx == Cont . (34) 


The last term in (383) represents the system performance degradation 
and is non-negative. 

The mean-square error during the early stage of data transmission 
can now be determined from (33) and (34). As the period of the training 
PN sequence approaches infinity, limy_. A;;(V) — 0, limy.. y.(V) — 0, 
and limy.. 6c — 0. Hence, C*, approaches C,,, asymptotically as 
N increases. 

We now assume some specific channel characteristics and apply 
these formulas to determine the initial performance degradation. 
Example 1: A baseband channel with a flat amplitude characteristic 
and a typical quadratic delay characteristic’ is assumed. The delay 
at the Nyquist frequency is taken to be 6,,7' seconds. The phase char- 
acteristic is of the form 


B(f) = 808,.T°f°/3. (35) 


The system impulse response, h(t), can be calculated from 
1/2T 
ni) = 2 | cos @nft + B()) af. (36) 
0 


In this example, we consider a typical value 8,, = 2. One hundred one 
samples of h(t) (from tj) — 507 to t, + 507) are taken with 7 = 1 and 
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t) = 0.67. For a 7-tap equalizer, the minimum mean-square error 
attainable is 0.03347. The mean-square error obtained from (83) is 
0.0393. The results for a 15-tap equalizer are’ 0.01484 and 0.01785, 
respectively. It is clear from these numbers that the performance 
degradation caused by using a PN sequence in the training period is 
negligible. This can be further illustrated by sketching the vector 5C 
in (84). Since the amplitude characteristic is constant, we have 

















A=I (87a) 
and 
A*=I-— 4, (87b) 
where 
1 1 1 
°N N N 
1 1 1 
Moy N 
1 1 1 
A = N N 0 N (37¢) 
1 1 
NNN : 
The inverse of A* is found to be 
. N-1 N? 
*)\71 
(A*) =I+tyayqityey: (38) 
From (10), (84), (37a), and (38), we have 
N-1 2N N’ 
ee =-yoivtyaitt wei i} t Ay. (39) 


Substituting (28), (32), and (87c) into (89), we obtain 
>» hi ww—1y 2) +t0+iN 


i~0 
8C = (40) 
De h_w-1)/21+to+iN 


740 


In Fig. 2, the time samples h, ,k = — © to ©, of A(t) are sketched. 
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ho 





Fig. 2—Illustration of the elements of 6C 


For N = 7 (7-tap equalizer), the first element of 6C is the sum of the 
infinite sequence --- , h_y, , h-4, Mio, «++ ; the second element of 8C 
is the sum of the infinite sequence --- , hy. , h_s , ho , «++ 3 ete. It 
can be seen that the large time samples h_; to hz are not included in 
these sums. This explains why 6C should have small elements. By 
repeating the above for N = 15, one can easily see that the performance 
degradation due to the use of PN training sequence is small and that 
this degradation approaches zero as N increases. 
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Example 2: We continue example 1 but change 8,, to its minimum 
value zero. The infinite sum ek h_xst.+in Can now be evaluated 
in closed form. By using the sampling theorem or the formulas for 
Psi (Digamma) functions,’ it is obtained after some manipulations 
ktN gin lo 


> hensto+in = (—1)" 


740 TW. 


N T a(ty — k) a(ty — k) 
{7 er ae (cis on. + 9 on )} oe 


where 0 S é S 0.57 is assumed. The performance degradation can 
now be determined in closed form 





N-1)]/2 32. 
Nee» gin tl 


2 2 
k=-(N-1]/2 aN 


{* E72 (cts sa + tg an m)} “) 


Smin and (8C)’A(SC) are plotted in Fig. 3 for N = 7, 15, and 31 and 
t) = 0.05, 0.1, 0.15, 0.2 and 0.25. It can be seen that the value 
of (6C)’A(6C) is approximately an order of magnitude less than that 
of Emin - Also note that (8C)’A(8C) reduces almost by half when N is 
doubled. These results again show that the performance degradation 
caused by using PN training sequence is negligible. 


(8C)’A(8C) = 





Vv. CONCLUSIONS AND DISCUSSIONS 


We have analyzed the start-up performance of a transversal equalizer 
for the case where a maximum-length pseudo-random sequence is 
used as the training signal to adjust the equalizer in the training period. 
The equalizer taps are adjusted by the gradient method to minimize 
the mean-square error, &*, between the equalizer output and the trans- 
mitted pseudo-random sequence. The pseudo-random sequence has 
been denoted 6; , B2, --* , By , 8: , Bo, °** , where N is the period of 
the sequence. We have considered the case where N is equal to the 
number of taps of the equalizer. The following results are obtained: 


(t) For a fixed N, the initial value of the mean-square error &*, 
the convergence rate of &*, and the minimum value of &* are 
all independent of the specific values of the 6,’s. Therefore, the 
same performance is obtained with any of the many pseudo- 
random sequences available. For example, a maximum-length 
pseudo-random sequence can be cyclic shifted to produce N 
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Fig. 3—Computational results of Example 2. 


pseudo-random sequences. Any one of these N sequences can 
be used in the training period with the same result. 

The initial value of &* depends on the phase characteristic of 
the communication channel, and the timing and phase settings 
at the receiver. However, the convergence rate of &* is inde- 
pendent of all these parameters. This result is similar to the 
one obtained previously* for the case where isolated test pulses 
are used as the training signal. 

Unlike the isolated test pulse case, the eigenvalues of the cor- 
relation matrix here are not close together. For example, when 
the channel has a flat amplitude characteristic, the first N — 1 
eigenvalues \, to Ay-; are equal to 1 + 1/N, while the Nth 
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eigenvalue Ay is equal to 1/N. The mean-square error can be 
decomposed into N components, each associated with one of 
the N eigenvalues. The tap gain adjustment reduces rapidly 
the N — 1 components associated with \, to Ay-1 , but the 
component associated with \y decreases very slowly. Therefore, 
as discussed at the end of Section III, care must be exercised 
in setting the timing, phase, and equalizer taps at the beginning 
of the training period so that the component associated with Ay 
has a small initial value. Note that this precaution is not re- 
quired when isolated test pulses are used, because in that case 
the eigenvalues are all close together and the components of 
the mean-square error all decrease rapidly.* 

The analysis shows that the tap settings obtained with 
maximum-length pseudo-random sequences with very short 
periods are nearly optimum for random data transmission. 
More specifically, the equalizer taps are adjusted in the training 
period to minimize the mean-square error &* between the 
equalizer output and the transmitted pseudo-random sequence. 
When such tap settings are used for actual data transmission, 
the mean-square error between the equalizer output and the 
transmitted random data can be written as 


& = &nin + € 


where Ein 18 the minimum attainable value of &, and e is non- 
negative, because tap settings obtained with pseudo-random 
sequence do not necessarily minimize &. Formulas for computing 
e were developed in Section IV and illustrated by numerical 
examples. It can be seen from these formulas and Figs. 2 and 3 
that « decreases rapidly as N increases (N is the period of the 
pseudo-random sequence and also the number of equalizer 
taps). The computations show that e decreases approximately 
by the factor 1/N (for example, « reduces approximately by 
half when N is doubled). The computations also show that « 
is about an order of magnitude less than &,;, . (This is so for 
N as small as seven.).Therefore, tap settings obtained with 
pseudo-random sequences are nearly optimum for actual data 
transmission. 

As can be seen from the computations in Section IV, &,in 1S 
rather large when N is small. For example, for a system with 
typical channel delay distortion (see example 1) and S/N = 
30 dB, Smin can be 15 dB above the thermal noise level when 
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N = 7, and 12 dB above the thermal noise level when N = 15. 
These large mean-square errors are due to the fact that for 
single-sideband Nyquist systems the overall system impulse 
response decays very slowly with time. Thus, for very sharp 
rolloff VSB systems (such as 4-percent rolloff), it is necessary 
to use a large number of equalizer taps (such as 31 or more). 


REFERENCES 

1. Lucky, R. W., and Rudin, H. R., ““An Automatic Equalizer for General-Purpose 
Communication Channels,”’ B.S.T.J., 46, No. 9 (November 1967), pp. 2179- 
2208. 

2. Mueller, K. H., Private Communication. 

3. Golomb, S. W., Digital Communications with Space Applications, Englewood 
Cliffs, N. J.: Prentice-Hall, Inc., 1964. 

4. Chang, R. W., “A New Equalizer Structure for Fast Start-Up Digital Communica- 
tion,’ B.S.T.J., 50, No. 6 (July-August 1971), pp. 1969-2014. 

5. Lucky, R. W., Salz, J., and Weldon, E. J., Jr., Principles of Data Communication, 
New York: McGraw-Hill, 1969, pp. 68-69. 

6. Gradshteyn, I. S., and Ryzhik, I. M., Table of Integrals Series and Products, New 


York and London: Academic Press, 1965. 


Copyright © 1972 American Telephone and Telegraph Company 
THe Bett SystTEM TECHNICAL JOURNAL 
Vol. 51, No. 9, November, 1972 
Printed in U.S.A. 


Stability Considerations in Nonlinear 
Feedback Structures as Applied to 
Active Networks 


By M. BAUMWOLSPINER 
(Manuscript received April 24, 1972) 


Active filters have recently acquired widespread use in the realization of 
frequency-selective networks. Unlike their passive counterparts, active 
filters have the potential of oscillating. 

Furthermore, it has been observed that the onset of oscillations in biquad 
active filters is dependent upon signal level. This led to the recognition 
that nonlinear stability theory would be necessary to comprehend this 
behavior. 

This paper develops a technique to analyze the stability of networks 
containing linear and nonlinear elements interconnected in multifeedback 
structures. This is accomplished by extending the concept of the ‘‘Describing 
Function’ to include networks containing nonlinearities with frequency- 
dependent linear feedback. The technique ts then applied to explain 
qualitatively and quantitatively nonlinear effects in op-amps and their 
relation to the stability of frequency-selective networks containing them 
(e.g., the Multiple Amplifier Biquad, MAB, and the Single Amplifier 
Biquad, SAB). The technique ts also applied to explain frequency shifts in 
amplitude-limited oscillators. The most valuable result of this analysis is 
the discovery of nonlinear feedback circuits which circumvent the conditional 
stability of high-frequency biquads. This has allowed us to obtain Q’s of 
50 at 100 kHz in a MAB employing 709 op amps. Similarly, a MAB 
employing 702 op amps was made to operate at 2 MHz with a Q of 10. 


I. INTRODUCTION 


In this paper, we shall deal primarily with a frequency-domain 
approach of analyzing networks containing linear and nonlinear elements 
interconnected as multifeedback structures. Particular applications 
will include the Single Amplifier Biquad’ (SAB), the Multiple Amplifier 
Biquad’ (MAB), and amplitude stabilized oscillators. 
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Fig. 1—A typical system with separable networks. 
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Our interest in this subject is motivated by conditional stability 
problems in biquads designed to operate above certain critical fre- 
quencies. Above these frequencies, these filters exhibit oscillatory modes 
when certain excitations are or were present. In the past, these 
conditional stability problems have been attributed to slew rate limiting 
in the operational amplifiers. However, we shall demonstrate later that 
although slew rate limiting may affect the stability, it is neither a 
necessary nor sufficient condition for conditional stability to occur 
in the MAB circuit. 

We shall draw heavily on the Describing Function (DF) and the 
Dual Input Describing Function (DIDF) Techniques.* ° These tech- 
niques have the virtue of imparting a conceptual understanding of the 
problem. Techniques based on the Liapunov stability criterion such as 
the Popov criteria and circle criteria have been looked into and seem to 
offer overly restrictive sufficient conditions for stability in most of our 
reasonably high Q applications. In addition, they fail to extend easily 
to multiple nonlinear loops. 

To facilitate the application of these techniques to the filter condi- 
tional stability problem, several interesting results concerning opera- 
tional amplifiers will first be derived. Earlier nonlinear op-amp models 
have represented an op-amp by saturating elements intermingled with 
frequency-dependent networks in a feedforward configuration. It will 
be seen that this is insufficient to predict insertion phase measurements 
of an op-amp in its nonlinear region. However, the usual position of 
compensation elements in an op-amp is in a feedback path around the 
gain stages yielding higher effective (Miller effect) capacitances. As 
a consequence, the linear and nonlinear elements become interconnected 
in feedback structures which impart varying phase characteristics as a 
function of amplitude. This becomes highly significant when considering 
active resonant circuits. Most important, when high Q circuits are 
realized, any additional phase lag around a loop may be sufficient to 
increase the pole Q to the point of oscillation. Similar phase shifts may 
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occur in oscillators in which nonlinearities are purposely introduced to 
stabilize their amplitude. In oscillators these effects are manifested 
by a discrepancy between the linearly-computed and actual frequency 
of oscillation. 

Finally, we will present circuits which produce either phase lag or 
lead as a function of amplitude. These circuits may be useful for non- 
linear compensation of the aforementioned problems. 


II. THE DESCRIBING FUNCTION TECHNIQUE 


The describing function method is an outgrowth of the Harmonic 
Balance technique used by Krylov and Bogoliubov,’ in nonlinear 
mechanics. The method reduces a nonlinear differential equation into 
a linear relationship by assuming a sinusoidal solution. The method is 
most useful when the system contains sufficient lowpass filtering to 
allow higher harmonics to be neglected. However, if the describing 
function technique is inadequate due to its neglect of higher order 
harmonics, the DIDF may help in solving the problem. 


2.1 Input-Output Concept of A Nonlinear Element 


We begin our analysis by reviewing some basic concepts of the 
DF technique. Consider the system of Fig. 1, where the linear and 
nonlinear parts are assumed separable. N is described in terms of its 
effect on a sinusoidal waveform. In particular, the describing function 
is defined as 


a Fundamentalofm M, 


Da) Fundamentalofe EF avis (1) 


In general, the DF will be a complex quantity. However, if N is a 
single-valued nonlinear function, its input-output characteristic will 
not enclose any area (see Fig. 2) and its DF will always be real as 
shown in Appendix A. 


m (QUTPUT) m (OUTPUT) 


e (INPUT) e (INPUT) 


(a) (b) 


Fig. 2—(a) A single-valued nonlinearity (real DF). (b) A nonlinearity having a 
iosteroais loop that encloses a given area (complex DF). 
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On the other hand, the DIDI is defined for a two sine wave input, 
with one wave at a multiple frequency of the first, i.e., 


= EF cos (wt + ¢) + EF, cos not. 


The DIDF is the ratio of the desired frequency component in the output 
over the same frequency component at the input. The DIDF has been 
worked out by West, ct al.,” for a saturating type of nonlinearity 
operating on a sinusoidal input in the presence of the third harmonic. 
This function can be used in determining the required perturbat’on 
of the DF when the third harmonic is not sufficiently low for the DF 
to be directly applicable. 

With the above ideas applied to Fig. 1, we may extend linear stability 
criteria to obtain the stable and unstable limit points for this case. 
Assuming that there is sufficient lowpass filtering following N in Fig. 1, 
the input-output relationship of the system is: 


vo _ _ 9(E)HGw) 


a Tae GHGs) @) 
_ Hi) a 
a + H(je) 


We now apply the Nyquist criterion, but instead of taking the critical 
point as —1 and incorporating g(/) into H(jw), we take —[1/g(E)] as 
the critical point. 

igure 3c shows the conventional Nyquist plot and the locus of 
—f[l /g({E)| with / as a parameter for the system the block diagram of 
which is given in Fi ig. 3a. The describing function of the nonlinearity,” 
as shown in Fig. 3b depicts the variation in amplitude (the phase shift 
being zero) of the fundamental of sin wt operated on by a dead-zone 
nonlinearity. We note that if # is small, the critical point is not encircled, 
and the system will be stable. As # is increased, we reach point A, and 
the system becomes unstable. As a consequence, F increases till we 
reach point B where the system enters into a stable limit cycle. The 
limit cycle is stable since if # increases, the system becomes stable 
causing I to decrease. On the other hand, if EH decreases, the system 
instability is such that F will increase. The intersection point B defines 
the amplitude (£) and frequency (w) of oscillation. It is to be noted 
that in the above case the frequency of oscillation occurs at the inter- 
section of H(jw) with the real axis since N introduces no phase shift 
in this case. However, when N introduces phase shift, this will not be 
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(c) 
Fig. 3—(a) System with dead-zone relay. (b) DF of dead-zone relay. (c) Nyquist plot. 


the case. Furthermore, even when N has a real DF and the second 
harmonic is not sufficiently filtered out, the DIDF predicts a phase 
shift of the fundamental through NV. 


2.2 Input-Output Concept of A Feedback Structure with Nonlinear 
Elements in tt . 


With the above as background, we proceed further with the DF 
concept by determining the sinusoidal input-output relationship for a 
feedback system with a nonlinearity. We will do this by considering the 
specific example of Fig. 4, which resembles an op-amp with unity 
feedback. As before, we assume sufficient lowpass filtering to eliminate 
the effect of the second harmonic. 
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Fig. 4a—System with saturating nonlinearity. 


It is clear from Fig. 4a that 
V;-V.=E or V,=E+YV,. (4) 


In addition, we know the relationship between EH and V, from the 
describing function of Fig. 4b and the linear transfer function H(S). 
Therefore, we can obtain the relationship between V; and V, . This is 
done most appropriately by considering the phasor diagram of Fig. 5b. 

First, we obtain the transfer characteristics of the linear block,. 
H (jw) as shown in Fig. 5a. Incidentally, this is the Nyquist plot for the 
linear region of operation. Assuming momentarily that we hold w = , , 
H(jw) and g(H’), determine the necessary magnitude and phase of H 
to give a particular V,. From eq. 4, we then find the corresponding 
V,; by taking the vector sum as is shown for several V, in Fig. 5b. 
Observe that the boundary between the linear and nonlinear region is 
(Vi?, VS, E™). In the linear region ZV, — ZV; is independent of 
amplitude, while in the nonlinear region ZV, — ZY, is a function of 
amplitude as shown in Fig. 5c. As a result of H(jw), there will be a 
different curve of this type for each frequency considered. These curves 
together with similar ones for | V, |/| V; | define the DF for the system 


aAe? an -9 cert [ane SIN ET ET (-et""| 






ZaiE'=0 


lE'| 


Fig. 4b—Describing function of saturating nonlinearity. 
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Fig. 55>—Vector diagram for determining V;. 
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Fig. 5c—Phase shift vs input signal. 


of Fig. 4a. This describing function differs slightly from those considered 
thus far in that it is a function of both amplitude and frequency. 

A most important result of the input-output concept developed in 
this section is the capability of being able to analyze a rather complicated 
nonlinear network by breaking it up into its constituents. Each one is 
analyzed individually and then combined as shown in Fig. 6. After 
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JEs[ SIN(w4 t+ /E,) » [Eg] SIN(w,t + /E2) |E3|SIN(w,t+ /E3) 
; : 


[Eq] SIN(w t+ /Ey) rd |Eg]SIN(w t+ /E3) 


93(E1,@7) = 9y(Eq,@7) 92[E191 (E1,@1),01] 


Fig. 6—Interconnection of DF’s 


reducing the network we apply a variation of the analysis employed 
for the example of Fig. 3. The difference lies in the fact that the DF is 
now a function of frequency and amplitude. As a consequence, we 
obtain many DF’s, each for a particular frequency, of which the negative 
reciprocal must not intersect the Nyquist plot at that frequency for 
stable operation. Usually, as we will see, the critical frequencies corre- 
spond to those close to the critical point in the Nyquist plane. 

A network which is well suited to this analysis technique is shown in 
Fig. 7. This is a typical Multiple Amplifier Biquad Filter using mono- 
lithic op-amps. Here, it is clear that if we employ the op-amp model 
shown in Fig. 7b, we have a multifeedback structure containing several 
nonlinearities. Deriving the DF of each closed loop op-amp and com- 
bining results, we can determine if the phase lag around the loop at 
certain amplitudes is sufficient to cause the biquad to oscillate. This 
is the subject of the next section. 


: OUT 
© 
— Aab 
( sta) (stb) 


Fig. 7—(a) MAB circuit. (b) A possible model for the op-amp. 
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Fig. 8—Fairchild 1A709 op-amp model. 


III. APPLICATIONS 


In this section we shall apply the techniques developed earlier to 
op-amps, active filters, and oscillators. Since active filters (e.g., MAB 
and SAB circuits) usually employ monolithic op-amps, we shall develop 
first an accurate model of the linear and nonlinear aspects of the op-amp. 
The op-amp model will generally depend on the type (e.g., 709, etc.) 
and compensation used. However, we shall demonstrate how to arrive 
at the model and present typical circuits. 


3.1 Operational Amplifiers 


3.1.1 Open Loop Characteristics of Op-Amps 


We have already presented one possible model of a typical op-amp 
in Fig. 7b. It is possible to modify the model by placing an amplitude 
limiter at the input or splitting the linear transfer function by inserting 
an amplitude limiter. An excellent model, as viewed from pulse and 
slew-rate measurements on a Fairchild nA 709 op-amp is shown in 
Fig. 8. The DF for this model does not predict any extra phase shift 
as a function of amplitude. This follows from Appendix A, where it is 
shown that this type nonlinearity (no hysteresis) introduces no addi- 
tional phase shift. However, if we consider the DIDI, we may indeed 
get extra phase lag as a function of amplitude. This results from the 
presence of harmonics, generated by the first nonlinearity, at the input 
of the second nonlinearity. The amount of phase lag may be estimated 





Fig. 9a—Symmetrical saturator. 


2038 THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1972 








Im (=— 





) 
DIDF 


Re (— 





DIDF 


a = AMPLITUDE OF 
FUNDAMENTAL 


b = AMPLITUDE OF 
THIRD HARMONIC 


RADIAL LINES ARE ALL 
10° APART AND 
REPRESENT LINES 
OF CONSTANT DIDF 
PHASE; ANNULAR LINES 
ARE LINES OF CONSTANT 
DIDF AMPLITUDE. 


Vo =V; -1<p; <1 
Vg =—~1 vi<-l 
Vo= 1 vi> 


Fig. 9o—Input-output relationships for fundamental in presence of third harmonic 


for saturating nonlinearity. 


STABILITY IN NETWORKS 2039 


by using the set of DIDF curves for a symmetrical amplitude limiter’ 
shown in Fig. 9. The curves for each value of a represent the DIDF as a 
function of the phase shift between the fundamental and third harmonic 
at the input of the saturator. 

T’rom these figures, it can be shown that a maximum fundamental 
phase lag of 11.1 degrees may be obtained when the op-amp of Fig. 8 
is heavily overdriven (i.e., a square wave input to the second saturator). 
Yet, when the phase characteristics of a 709 op-amp were measured 
in the lab, no such effect was observed. The curves obtained are shown 
in Fig. 10. These curves display phase lead as a function of amplitude 
instead of phase lag. This result is extremely important, for it verifies the 
stability inherent in biquads at low frequencies as we shall see later. 
In summary, a better model of the op-amp is needed. 

To determine this model, we first consider a simple circuit depicted 
in Fig. 11. This circuit is a transistorized amplifier which provides 
output limiting of the signal. Its applications include FM limiters and 
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Fig. 10—Measured phase characteristic of an open loop 709 op-amp. 
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Fig. 11b—Model of Fig. 11a assuming transistor is in linear region. 


Ig(Vo) 





Fig. 11c—Characteristics of diode pair. 


amplitude-stabilized oscillators. Defining Z as 


= — fe __ 
~ SCR, +1 


it follows from the circuit that: 


[= V/Z + Ti(Vo) (6) 


Z (5) 


STABILITY IN NETWORKS 2041 





Fig. 12a—Block diagram of Fig. 11b. 


or 
Vo = 2 — I.(Vo)]. (7) 


This equation is represented by the block diagram of Fig. 12a. Knowing 
the relationship between J, and J;, we proceed to draw the phasor 
diagram of Fig. 12b. Initially at low levels, 7, is small compared to J, 
because the diodes are cut off. However, if Z, is just short of turning 
the diodes on, a slight increase in J, will cause a large increase in J, . 
These two situations are represented by (Z§”, If?) and (I, I) 
respectively, in Fig. 12b. Taking their vector sum, we obtain J“ and 
I. Interestingly we note that the phase difference between J and I, 
and hence V;, and V, has decreased. By taking the small signal case as 
a reference, we conclude that the closed loop network exhibits phase 
lead with increasing amplitude. 

A significant conclusion to be drawn from this example is that although 
the simple circuit appeared to contain just a clipper (which would 
indicate no amplitude-dependent phase shift) the feedback present 
alters the situation sufficiently to predict phase shift as a function of 
amplitude. 

If the same form of analysis is carried out for the internal circuitry 
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Fig. 12b—Input/output relationship. 


2042 THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1972 


, SLOPE = —p=~9,.R3 







SC,R> 





SC, (R,+R5+R,)+1 








0.7 VOLT 


(a) 






R3 





SC, (R,+ Ry) +1 





SC, 






(b) “S SLOPE k 


Fig. 13—709 op-amp input compensation circuity model. (a) Cutoff region. (b) 
Saturation region. 


of the 709 op-amp, similar results are obtained. This has been done in 
Appendix B with particular attention to the circuitry involved in the 
input and output compensation. Most often, compensation in an op-amp 
is obtained by making use of the Miller effect to obtain low-frequency 
breakpoints with relatively small capacitors. Both the input and output 
compensation in a 709 op-amp take advantage of this effect. As a 
result, the transistor nonlinearities have RC feedback around them. 
The equivalent circuits of the input compensation circuitry in the 
cutoff and saturation regions, as derived in Appendix B, are shown in 
Tig. 18. In this figure we have left out some linear transfer function 
blocks at the input and output side, since these have no effect on the 
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Fig. 14a—DF for nonlinear element in cutoff region model. 
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nonlinear characteristics. The following analysis will show that both 
of the circuits in Fig. 13 provide phase lead as amplitude increases. 

Figure 14a shows the DF of the nonlinear element’ in the cutoff 
region model under the assumption that the quiescent point is at 
E = EH, . Figure 14b depicts the phase relationships of the linear element 
in that model. With this information and Fig. 13a, we construct the 
phasor diagram of Fig. 14c. Specifically, since the nonlinear element 
is single valued and therefore contributes no amplitude sensitive phase 
shift, —V, will lead H according to the linear element and the inversion 
due to the nonlinear element. Entering the nonlinear region of the 
model, the error voltage £ has to increase at a faster rate than V, to 
overcome the attenuation effect of the clipping element. From the 
relationship 


Vie (Vg) ek (8) 
AVo— %V, 
90° 
45° 
@ 
C,(R,+R5+R3) 


Fig. 14b—Phase characteristic for linear element in cutoff region model. 
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Fig. 14e—Phasor diagram for the closed loop in cutoff region model. 
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Xo = Voo-Vsat. 


Fig. 14d—DF for nonlinear element in saturation region model. 


we determine the third vector V,. As a result, the angles between 
(—V5P, Vi?) and (—Vi”, V%) do change with increasing drive. 
However, at high frequencies where H and —V>, are initially close to 
being in phase, the increase of phase lead of V, relative to V; is minimal. 
The frequency breakpoint below which the nonlinear phase lead effect 
may be expected is typically, 


fa ee 
oa 2n(hy + Rh. + R3)Ci 


On the other hand, in the saturation region, the model of Fig. 13b 
applies. Here the DI for the nonlinear element is given in Fig. 14d under 
the assumption that the quiescent point is at Vo = Voo. Furthermore, 
we have approximated the exponential nonlinearity of the diode in a 
piece-wise linear manner. It follows directly that the phase relationship 
of H with respect to V, is given by the phase characteristic of the lag 
network following the nonlinear element in Figure 13. Drawing the 
phasor relationships for the closed loop at a fixed frequency where the 
lag network has singificant phase lag yields Fig. 14f. Again, the phasor 
diagram reveals that the phase lead of V, with respect to V; increases 
as amplitude increases. The frequencies over which this phase lead 


~ 50 kHz. (9) 
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Fig. 14e—Phase characteristics for the linear element in saturation region model. 
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Fig. 14f—Phasor diagram for the closed loop in the saturation region model. 


effect will be present is approximately the range of frequencies where 
the linear network provides phase lag. From Fig. 14e, this occurs for 
frequencies below: 


2 


i! 2nC,(R, + Ry) ~ 100 kHz. (10) 


Above this frequency range, #‘” and EF will have approximately the 
same orientation as V{? and V{”, forcing V{? and Vi” to maintain 
their orientation roughly independent of amplitude. 

The above analysis also applies to the output compensation circuitry 
provided that it is properly interpreted. The output compensation cir- 
cuitry for a 709 op-amp and its nonlinear model are shown in Fig. 15. 
Here, unlike the input compensation, the RC network is fed back from 
an emitter follower circuit which simplifies the model and increases the 





Fig. 15a—Output compensation of 709 op-amp. 


2046 THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1972 


_-SLOPE =— 





gts) ee 
S+ 1/C, Rp 


Fig. 15b—Nonlinear model of 709 op-amp. 


breakpoint frequencies. In addition, the RC network is fed back to a 
lower impedance termination which also has the effect of increasing 
the breakpoint frequencies. In a typical 709 op-amp, this breakpoint 
frequency occurs at (the equivalent of eq. 9 with R, and R; being zero) 
pigs 25 8S 6 
{re nC R, © 1.5 X 10° Hz. 
This frequency is significantly higher than those resulting from the 
input compensation (see eqs. 9 and 10). As a result, the output compen- 
sation will be quite irrelevant to the conditional stability of biquad 
filters in the 0 to 100 kHz range. Above 100 kHz, the amplifier will no 
longer be used in a biquad filter, since the open loop gain is less than 
40 dB, which is insufficient for most precision applications. 

We now have an op-amp model consistent with the experimental 
data of Fig. 10. Having discussed these open loop characteristics of an 
op-amp, we shall next use these concepts to derive the phase properties 
associated with closed-loop operational amplifiers. 


3.1.2 Closed-Loop Characteristics of Op-Amps 

We shall confine ourselves, in this section, to gain inverters, inte- 
grators, and leaky integrators. The circuit of Fig. 16a is the general case 
which includes all the above-mentioned configurations. The model of 
the closed loop op-amp shown in Fig. 16c is derived as follows. Assuming 
that the op-amp has infinite input impedance and letting 

€in = O, 
the voltage developed across the input of the op-amp will be 
R, SRC, + B,/R, 
eé= 1 out = SRC, “ [R,/R, zs 1] out + (11) 
R, + Ea || R, | 
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Fig. 16a—Closed loop op-amp. 
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Fig. 16c—Model for Fig. 16a. 


Likewise, letting 
Cout > 0, 


the input voltage of the op-amp will be 


1 
SC, || Ra _ 1 : 
SRC, + [R,/R,g + 1" 





¢= (12) 


Se ig 
i 
Ry - Ea | R, | 
The block diagram of the closed-loop op-amp, obtained through the 
superposition of eqs. 11 and 12, is shown in Fig. 16b. By the reduction 
method, we move the feedback network past the summing node and 
make the proper correction to the input network. This leaves us with 


the model of Fig. 16c. A significant advantage of this model is that 
the ideal transfer function has been separately realized by the network 
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between V/ and V;. Consequently, the feedback network following 
the ideal transfer function network gives a direct measure of the error 
introduced by the nonideal characteristics of the op-amp. 

For the model of the op-amp itself, we may either use the networks 
of Fig. 13 in a cascade combination, or we may directly employ the 
DF curves shown in Fig. 10. The latter method is to be preferred since 
it is based on actual measurements and requires a minimum of computa- 
tion. Fig. 17 is the phasor diagram for an inverter (i.e., Fig. 16 with 
7 = Oanda = 1) at f = 5 kHz. In the linear region, we initially obtain 
a very small error voltage H and consequently V; and V, are almost 
equal. This error voltage, although small, has a distinct phase lead of 
approximately 70 degrees with respect to the output voltage V,. 
This phase shift is introcuded by the frequency characteristics of the 
op-amp and the network N,(s) in Fig. 16c. However, for the unity gain 
inverter, N.(jw) is a real quantity and, therefore, will not add additional 
phase. As we enter the nonlinear region, the magnitude and phase of 
the error voltage, H change as dictated by the describing function (DF) 
of the op-amp. The phase changes according to Fig. 10 and the amplitude 
according to the real DF of a saturation type of nonlinearity shown in 
Fig. 4b. Although the op-amp does not behave exactly as a saturation 
nonlinearity; the magnitude of the DF for the op-amp is well approxi- 
mated by this DF, unlike the phase shift of the DF which for the satura- 
tion nonlinearity is zero. 

Having obtained the phasor diagram of Figure 17, we note that 
increasing the input V; predicts an increasing phase lag of the output 
—V, relative to the input V,; . This effect has been verified experimentally 
as shown by the curves in Fig. 18. At 27 kHz where the op-amp has 
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Fig. 17—Phasor diagram for unity gain inverter (using a 709 at 5kHz). 
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Fig. 18—Measured phase characteristics of a unity gain inverter employing a 709 
op-amp. 


reduced gain, the effective saturation nonlinearity comes in earlier 
producing the phase lag at a lower input level. 

For the integrator, the same analysis applies but with the critical 
difference that N2(jw) in Fig. 16c is no longer a real quantity as in the 
case of the inverter. The phasor diagram for this case is shown in 
Fig. 19. We have chosen an integrator with a = 0 and 7, the reciprocal 
of angular frequency, 27-5 kHz. As a result, N2(jw) will provide a 
45-degree lead at 5 kHz, which will have the effect of moving the phase 
of the error voltage, H, 45 degrees closer to the vector V, than in the 
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Fig. 19—Phasor diagram for integrator at 5kHz, 7 = 1/275kHz, a = 0. 
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case of the inverter. However, when the system is overdriven, / swings 
around to the lagging side of V, , which consequently causes V, to lead 
V,. Asa result, in the integrator at this frequency, increasing the drive 
level produces a phase lead effect. This has been verified experimentally 
and the results are shown in Fig. 20. 


3.2 Active Filters 


3.2.1 The Multiple Amplifier Biquad Filter 


The above discussion has provided a basic feel for what happens when 
the op-amp is employed in a simple feedback loop. We shall next discuss 
the performance of networks which use these inverters or integrators 
as building blocks. 

The MAB’ circuit was already shown in Fig. 7a in its bandpass 
configuration. The configuration generates a pair of poles with arbi- 
trary Q and frequency location. To generate a pair of zeroes, outputs 
taken at different points are summed in a separate summing amplifier. 
The zeroes can also be generated by feeding the input to more than 
one op-amp as is done in the Multiple Input Biquad Filter (MIB). 
However, their basic operational frequency limitation manifests itself 
by a conditional stability problem in its pole forming loop. This problem 
is more acute when high Q’s are desired, for then the resonator loop 
operates very close to the critical point of its Nyquist plot. This is 
shown in Fig. 21, where it is seen that the phase shift around the loop 
at the pole frequency gets very close to 180 degrees as the Q is increased. 

As a consequence, a small amount of additional phase lag in any of the 
amplifier blocks (i.e., the inverter and the leaky and nonleaky inte- 
grators) making up the main loop can make the system unstable. 

With the information already acquired, as to the characteristics of 
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Fig. 20—-Measured phase characteristics for integrator employing a 709 op-amp 
(the gain is unity). 
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Fig. 21—Nyquist plot for main loop in MAB. 


the individual amplifier blocks, the task of determining whether a 
given filter will be stable can be investigated. By adding up the phase 
shift curves of the individual amplifiers, as determined from Figs. 18 
and 20, we are able to determine if the outer loop is stable with increasing 
drive level. As an example, we shall consider a MAB circuit which has 
a pole frequency at 5 kHz with a Q of 20. The design is such that 


1 i 1 


pe RO 0 
- WR;R2C,C; RC, RC, " 
and. 
_ i 
Q = re 20 


In this case, the integrators will have approximately unity gain at the 
resonant frequency of the pole. Consequently, all the three amplifiers 
enter the nonlinear region at roughly the same level. From the given Q 
we determine the phase lag needed to bring the main loop into oscillation 
at the pole frequency. For the reasonably high Q of 20, this is given by: 
on = a 0.05 rad = 2.9°. 
Q 

Referring back to the phase characteristics of the inverter in Fig. 18, 
we note that in the inverter, phase lag is introduced as the input level 
is increased. However, this is counterbalanced by the phase lead intro- 
duced by the integrator (see Fig. 20). 

It is to be noted that for very high input levels where the phase lag 
of the inverter becomes dominant, the circuit will, nevertheless, be 
stable. This results from the sharp drop in the magnitude of the DF 
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when the amplifiers are heavily saturated. The drop in the magnitude 
of the DI has the effect of shifting the critical point toward the left 
in the Nyquist plot of Fig. 21. 

If the same circuit is scaled up to 27 kHz, the outer loop will un- 
doubtedly be conditionally stable. In this case, both the inverter and 
integrator provide phase lag (Figs. 18 and 20) and at the level which 
produces 2.9-degree phase lag, the circuit will become unstable. As the 
signal level increases, due to the instability, the circuit will reach a 
stable limit cycle where the magnitude of the DIF has decreased suffi- 
ciently to just touch the Nyquist curve. 

The method of analysis employed thus far has its greatest advantage 
in that it is able to suggest many ways of getting around the stability 
problem in MAB filters. These methods center in either improving the 
manner in which the magnitude of the DF drops (possibly by making 
this occur earlier so that it becomes dominant) or improving the phase 
characteristics of the DI’. One method which has been employed 
successfully is the use of a diode network across the leaky integrator. 
This has the effect of improving both of the DF characteristics, magni- 
tude and phase. This circuit and the model for it are shown in Fig. 22. 
The model in Fig. 22c is derived very simply by taking 


t = v and 0 volt at op-amp input. 
We then have, 


Cin 





i= er = ise) + eoSC, + 5 =i! 









igleo) DIODE 
iz NETWORK x 





(a) (b) (c) 


Fig. 22—(a) Diode network across integrator. (b) I/O characteristics of the diode 
network. (c) Model for (a). 
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which if rearranged, we get the equation describing the model 


—€in ' 1 
€y) = | Ro isle) | = 
: SC, + R, 
Q 





Drawing a phasor diagram for this circuit reveals that the magnitude 
of the DF drops and simultaneously phase lead is introduced as signal 
level increases. It is worthwhile to point out that limiting directly at 
the output of the op-amp has the opposite effect (as shown earlier in 
the example of Fig. 4); therefore, limiting cannot be blindly applied. 
This technique has been employed in a MAB using 709’s to obtain 
a Q of 50 at 100 kHz. Above this frequency, the gain of a 709 op-amp 
is no longer sufficient to give a precision filter. Similarly, a MAB em- 
ploying 702 op-amps was made to operate at 2 MHz with a Q of 10. 

This method of nonlinear compensation has the disadvantage of 
limiting the dynamic range of the filter. Therefore, techniques which 
compensate the op-amps internally are more desirable. It has been 
shown in Section II how internal compensation. affects the DI of 
closed loop amplifiers. In turn, these circuits may be rearranged to 
produce better characteristics. We will not dwell on this subject in great 
detail; however, we shall discuss the connection with slew rate. It is 
well known that slew rate limiting in operational amplifiers is caused 
by some form of either voltage or current limiting. Consequently, as 
a result of this limiting, we may expect that the DF and, therefore, 
conditional stability will be affected. However, this may not necessarily 
be the case for it may happen that the limiting clement causing slew 
rate is not the dominant factor in determining the phase characteristic 
of the op-amp. This is most vividly illustrated by a 709 op-amp with 
input and output compensation. It has been experimentally observed 
that the input compensation affects the frequency range of conditional 
stability while the output compensation has a negligible effect. On the 
other hand, slew rate limiting is mostly affected by the output compen- 
sation. The dominance of the input compensation in controlling condi- 
tional stability is consistent with our derivation of the DF of an open 
loop 709 op-amp.” 


3.2.2 The Single Amplifier Biquad Filter 


The single amplifier biquad’ circuit is shown in Fig. 23a. In Figs. 23b 
and c, we have the feedback structure, as seen by the op-amp, and its 
Nyquist plot. Here, once again, we deal with a single op-amp in a closed 
loop configuration and our interest lies in determining its stability. 
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Fig. 23—(a) SAB Circuit. (b) Feedback structure as seen by op-amp. (c) Nyquist 
plot of feedback structure assuming Rj|| Re > Re. 


From the open loop DI curves of the op-amp, we can plot the locus 
of the critical points as —1/DF. As shown in Fig. 28c, the circuit is 
conditionally stable as a result of more than 90 degrees of phase lag 
in the amplifier. One manner in which stability can be achieved is by 
appropriately including a diode network in the feedback structure such 
as to move the circle in Fig. 23c to the right when the drive level exceeds 
a predetermined threshold. This network can cither be placed across 
R, or Rp . Another alternative, which is advantageous from a dynamic- 
range perspective, is to design the op-amp or its compensation such 
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that proper phase lead is introduced sharply as the level is increased. 
Modified forms of the compensation circuitry discussed in Section 3.1 
can be used to achieve this effect. 


3.3 Amplitude Stabilized Oscillators 


Sine wave oscillators are usually designed in the linear domain by 
placing the poles slightly to the right of the jw-axis in the S-plane. If 
the network remained in its linear region, the network output would 
increase exponentially without bound. Invariably, therefore, nonlineari- 
ties are introduced to produce a stable limit cycle. Fig. 24 illustrates 
a typical circuit which makes use of a saturation type of nonlinearity. 
Since the DI of a saturating nonlinearity is real, the circuit will os- 
cillate at the frequency of the tuned circuit, independent of the non- 
linearity. 

Many times, especially at high frequencies, nonlinearities with an 
effective real DF are difficult to obtain. In these cases, the computation 
of the DF (as a function of amplitude and frequency) is required to 
determine both the amplitude and frequency of the limit cycle. To 
illustrate this point, we refer to the oscillator of Fig. 25. (This circuit 
is being currently employed in a 20-MHz subcable oscillator. In this 
application, reliability necessitates an accurate determination of the 
frequency of oscillation.) A first look at the circuit may lead us to the 
conclusion that it behaves in the same manner as the oscillator shown 
in Fig. 24a. The pair of diodes would behave as the saturating non- 
linearity while the frequency selectivity would be provided by the 
tank circuit. However, a second look reveals that this circuit is quite 
similar to that shown in Fig. 11. It can also be seen, recalling Fig. 12, 





(a) (b) 
Fig. 24—-(a) Amplitude-stabilized oscillator. (b) Nyquist plot. 
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Fig. 25—Amplitude-stabilized oscillator by the use of diodes. 


that the pair of diodes is effectively in a feedback configuration with 
frequency-sensitive elements, which produces a complex DF. Therefore, 
we conclude that this circuit will have an oscillation frequency depen- 
dent on the nonlinearity. 


IV. CONCLUSION 


A technique which analyzes the stability of frequency-selective 
feedback structures containing nonlinearities has been presented. This 
technique is an adaptation of the describing function technique to 
include multiple feedback structures. This is accomplished by extending 
the definition of the conventional describing function to include net- 
works containing nonlinearities around frequency-dependent linear 
elements. By doing so, the describing function of typical operational 
amplifiers in open and closed loop configurations have been derived. 
These, then, serve as a means of predicting the conditional stability 
criteria in frequency-selective feedback networks such as the Multiple 
Amplifier Biquad and the Single Amplifier Biquad. A valuable result 
of this analysis is the discovery of circuits which circumvent in one way 
or another the conditional stability of high-frequency biquads. This 
result has had a highly beneficial impact on high-frequency biquads 
since it relieves the conditional stability problem and shifts emphasis 
to the maximum open-loop gain available at high frequencies. This 
open-loop gain determines the precision which high Q circuits can 
achieve at these high frequencies. 
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APPENDIX A 


A Property of the Imaginary Part of A Describing Function 


Theorem: The imaginary part of the describing function, G(A), associated 
with a multiple valued input-output function, f(a), is given by 


=e 

wA® 

where S, ~ area enclosed by f(a) for |a| < A; positive when path taken 
by nonlinearity 1s in a counterclockwise direction. 


Im {G(A)} = 


Proof: 
Im {@(A)} af f(A sin wt) cos wt dw). (13) 
Let, 
= A sin at. 
Then, 


dp = A cos ot d(wt). 
It follows by substitution into eq. 13 that, 


qe? 1Wdu= -St QED. 








Im {G(A)} 


APPENDIX B 


709 Op-Amp Models 


The schematic of the input circuitry of a 4709 op-amp is shown 
in Fig. 26. The transistors Q, and Q, make up, the second stage of 
amplification. Q;, Q;, and Q, decode the differential output of Q, 
and Q. providing a single-ended signal at the base of Q,. The input 
compensation, a series RC network, is placed across base and collector 
of the Darlington pair, Q, and Q,. This provides (together with the 
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Fig. 26—Input circuitry of uA 709 op-amp. 


output compensation) the necessary frequency rolloff to stabilize the 
amplifier in the linear domain. 

Fig. 27 depicts the approximate circuit model of the second-stage 
amplifier derived from the Ebers—Moll model.*’ We have assumed that 
the input impedance of the Darlington stage is much larger than R, 
and the load impedance presented by the emitter follower, Q,;, is 
negligible compared to R;. The cutoff region is represented by the 
characteristics of V,, shown in Fig. 27b. Saturation is introduced in 
the model by the diode D. In Fig. 27a the current I is given by 


Cin oe Vo Vo VCVe.) 


I = Se ot Va) 
R, + Ri + 


A 
SC, 
Rewriting this equation, gives 


Z Zo 
Vo = & ein + Val Vee) B — to(Vo)Zo (14) 
Zi, Rs 
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where, 









SLOPE = — u=-9,,R 
Vp 


e 


Fig. 27b—Nonlinear elements of the second stage amplifier corresponding to the 
cutoff region. 





Fig. 27c—Nonlinear elements of the second stage amplifier corresponding to the 
saturation region. 
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We also obtain from Fig. 27a 


Vie = Cin — (Coan, 


Z, 
= Cin Ry = Vo Es (15) 
3 4“, 
where, 
_ 4iho 
ae eae 


Figure 28 is a block-diagram representation of eqs. 14 and 15. For 
simplicity, we distinguish between the cutoff and saturation regions. 
The result of this is shown in Fig. 29. In both models, the forward path 
containing Z,/Z, in Fig. 28 is neglected since the other path through 
V, has a much larger gain. We have also assumed, with some loss of 
generality, that the biasing is such that a sufficiently high drive level 
producing cutoff (saturation) does not cause saturation (cutoff) during 
some other portion of the cycle. However, in general, the problem can 
be solved by computing the DF for the inner loops involving 7p and V, 
and employing these results to solve for the DF for the complete 
system. 

The model of Fig. 29a represents the second-stage arnplifier when 
driven into cutoff. It is obtained directly from Fig. 28 by making a 
slight rearrangement at the output in addition to the two previously 





Fig. 28—Block diagram representation of eqs. 14 and 15. 
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Fig. 29—(a) Model for the cutoff region. (b) Model for the saturation region. 


mentioned assumptions. The network parameters are given by 


ZR, SCR, 


ZR;  SC\(R, +R. + Rs) +1 


Ry SCR, +1 
Z, SCR, +R.) +1 


Z, _ SCR + Ry) +1 
Re SCR, 


The model shown in Fig. 29a depicts the saturation region. It is obtained 
from Fig. 28 by restricting V, vs V,. to operate in the linear region. 
Consequently, this nonlinear element (V, vs V,.) in Fig. 28 becomes a 
linear gain element. The process of reduction is shown in Fig. 30. Figure 
30c, after factoring out the feed-forward path and assuming a large up, 
yields the saturation model of Fig. 29b. The linear element in the feed- 
back loop in terms of fundamental parameters is: 


RZ, _ Rs | Sete + R:) + ih 
Roy = Rop SC, 
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Fig. 30—Reduction of Fig. 28. 
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Waveguide Breakdown Effects at High 
Average Power and Long Pulse Length 


By A. S. ACAMPORA and P. T. SPROUL 
(Manuscript received March 17, 1972) 


Analytical predictions of the power handling capabilities of waveguide 
systems generally have not considered the effects of high average power and 
long pulse length. It has experimentally been noted, however, that a sub- 
stantial reduction in power handling capability below expected levels does 
occur as average power and pulse length are increased. This reduction can 
be attributed to the presence of loose particulate matter which is heated 
by the average power, causing localized rarefication of the dielectric gas 
fill as well as the expected voltage enhancement. 

An unstable arcing situation is shown to exist when the arc duration 
exceeds some critical time. Typical pulse lengths in use today exceed this 
critical time and may result in continuous arcing. The use of control 
circuitry to terminate each particle-induced arc prevents continuous arcing 
and deletes the particle, and is therefore essential to stable operation at 
long pulse lengths and high average power. 


I. INTRODUCTION 


An understanding of microwave breakdown in gases has existed for 
many years’’” and has been successfully applied to the design of many 
high-power microwave systems operating at low duty cycles with pulse 
lengths on the order of several microseconds. However, recent testing 
of S-band WR284 microwave hardware filled with sulfur hexafluoride 
(SF,) and operated at 500- to 1500-kilowatt average power, and 10- to 
150-megawatt peak power with 50- to 150-microsecond rectangular 
pulses showed performance far below expectation. A series of tests and 
supporting analytical work were undertaken to explain this behavior. 
These have shown that the empirically observed reduction in peak power 
handling capability at high levels of average power is caused by the 
presence of loose particulate matter, microscopic in size and lossy at 
microwave frequencies. The experimental and analytical proof of the 
mechanisms responsible for this effect is described. A procedure is 
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presented employing short duration discharges to break up these 
particles, thereby increasing power handling capability. 

Waveguide cavities and resonant ring structures of various lengths 
were tested because the empirical peak breakdown power level was 
found to decrease with resonator length. It is shown that this effect is 
related to the arc energy absorption which is dependent upon resonator 
length. A quantitative understanding of this phenomenon is developed 
and is applied to confirm the experimental results. 


II. BREAKDOWN MECHANISMS 


The parameters which affect the peak power handling capability of 

a microwave system are: 
(¢) the maximum electric field strength appearing within the 
system 
(21) the nature of the dielectric gas employed 
(121) the molecular density (or, equivalently, the absolute pressure 
and temperature) of the dielectric gas 
(iv) the microwave frequency 
(v) the breakdown volume. 
Free electrons which always exist within the gas because of cosmic 
radiation and other random phenomena are accelerated by the electric 
field present and suffer collisions with neutral gas molecules. If the 
kinetic energy imparted to an electron by the field exceeds the ionization 
potential of the gas, the possibility exists of producing an additional 
electron upon collision with a neutral gas molecule. If the kinetic 
energy of the accelerated electron is less than the ionization potential, 
it may become attached to a neutral gas molecule upon collision. 
Diffusion produces a net flow of electrons from regions of high electronic 
density to regions of lower density. Some recombination of electrons 
and positive ions also occurs. 

Breakdown occurs when the rate of electron production via ionization 
exceeds the combined rate of electron loss through the processes of 
attachment, diffusion, and recombination. lor most high power trans- 
mission systems, the gas pressure is in the region of one to two atmos- 
pheres, and recombination rates are negligibly small. Also, since the 
dimensions of the breakdown volume are generally much greater than 
the characteristic diffusion length of the gas, electron loss via diffusion 
is also negligible. Hence, breakdown occurs when the rate of ionization 
exceeds the rate of attachment. 

Both ionization and attachment rates are functions of the mcan 
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electron energy® which, in the high pressure regime where collision 
frequency is far larger than the microwave frequency, is dependent on 
the ratio of the peak electric field strength (/) to the gas molecular 
density (NV). Assuming the density is defined by the ideal gas law, the 
field strength at breakdown is given by 


E=CER (1) 


where 


p = absolute gas pressure 
IX = Boltzman’s constant 
T = absolute temperature 
C = constant, dependent upon dielectric gas 


Tor SF, , C = 3.69 X 107” volt-meter*/molecule. 


For irregularities in the waveguide such as those which occur at 
flanges, in hybrids, etc., the local voltage gradient is increased by a 
field enhancement factor (@). Substituting into (1) the waveguide 
field strength required for breakdown becomes 


1 = BKT (2) 
The electric field strength generated in a rectangular waveguide 
transmitting a power P can be shown to be 


_ 4Z,P 


2 
i= qe. (3) 





where d, and d, are the waveguide cross-section dimensions and Z, is 
the wave impedance. Therefore, the breakdown power threshold can 
be derived from (2) and (8) as 


did C” 2 
si (482)(sora)e (4) 


Tor constant temperature, this is the familiar relationship between 
peak power breakdown and the square of the absolute pressure, known 
to be valid in the high-pressure region where electron loss by diffusion 
is negligible. 





Ill. EFFECTS OF HIGH AVERAGE POWER 


The relationship discussed above has been understood and effectively 
applied for many years. However, in the experimental work described 
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below, a significant reduction in peak breakdown power was observed 
as the average power increased. It is noted that, although eq. (4) is not 
explicitly dependent upon average power, temperature does appear. 
Hence, it was postulated that loose particulate matter was reducing 
the breakdown power level because a loose particle is not thermally 
attached to the waveguide wall; its temperature therefore increases 
as the average power increases. This causes a reduction in the molecular 
density of the gas at the location of the particle, thereby reducing the 
field strength needed to cause breakdown. The presence of a particle 
also produces localized field enhancement. These two localized effects 
combine to lower the breakdown power. 

The electric field enhancement factor and the equations for tempera- 
ture rise at various levels of average power are derived in Appendix A 
for spherical conducting particles. This derivation neglects intrapulse 
heating since, within the range of pulse lengths investigated, it is 
negligible compared to average power heating. Using these results 
and eq. (4), peak power breakdown levels can be calculated for various 
average power levels and particle sizes. Figure 1 shows the results of 
such calculations for WR-284 waveguide filled with SI, at 25 psig 
containing spherical copper particles from 4 to 120 microns in diameter. 
Also appearing on Tig. 1 are experimental points obtained by intro- 
ducing copper spheres of known diameter into a traveling wave reso- 
nator’ constructed of high-conductivity copper to minimize losses. 
Correlation between theory and experiment is seen to be good. 

It is observed from Fig. 1 that, within the range of particle sizes 
investigated: 


(z) Peak breakdown power levels at very low levels of average power 
are independent of particle size and are below that of ideal 
waveguide because of field enhancement only. 

(27) As the average power increases, the presence of heated particles 
causes further reduction in power handling capability due to 
localized reduction in gas density. Larger particles cause a 
greater reduction in power handling capability because they are 
heated to higher temperatures. 


The analysis from which Fig. 1 was derived is applicable when the 
particle dimensions are small compared with the waveguide dimensions, 
and large compared with the characteristic diffusion length of the gas. It 
should not be inferred from Fig. 1 that particle-induced arcing cannot 
occur above 100 megawatts, since particles with dimensions smaller 
than the diffusion length will cause breakdown to occur above this level. 
An analysis of particle behavior in this region was not considered. 
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Fig. 1—Effect of particles on breakdown. 


IV. EXPERIMENTAL WORK 


The equipment used to obtain the experimental data is shown in 
schematic form in Fig. 2. A resonant ring was used to generate the 
high power levels required to produce breakdown with the RF power 
source available. A photograph of the equipment is shown in Fig. 3. 
View ports were provided to permit visual observation of arcing. 
Thermocouples were inserted in flanges and attached to the waveguide 
in various places to monitor temperature rise versus average power. 
Forward and reverse power was monitored by power meters and by 
oscilloscopes. A counter was used to record the number of arcs. 

Because of the high average power levels, the waveguide was water- 
cooled. Two cooling channels were soldered on the top and bottom walls 
through which water was circulated at a rate of 7 gallons per minute. 
To permit assembly of the flange bolts, the cooling channels were 
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Fig. 2—Resonant ring experimental setup. 
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Fig. 3—Resonant ring. 


terminated 2-1/2 inches from the flange. Because of this, the flange 
temperature was found to rise at a rate of 65°C/megawatt of average 
power. This factor has been included in the particle temperature calcu- 
lations of Appendix A. The flange used was adapted from a vacuum- 
tight design, the principle feature of which was the use of a stainless 
steel knife edge and a soft copper gasket to provide the RF and gas seal. 
Figure 4 shows a model of the flange and cooled waveguide. 


V. RF PROCESSING 


Early experience with the setup described above showed arcing at 
substantially lower peak power levels than expected. By use of the view- 
ports, it was observed that many arcs were being initiated at the power 
coupler and monitoring coupler slots. These were attributed to high 
gradients at sharp corners. By allowing arcing to continue, the test 
setup was “aged” to progressively higher and higher peak and average 
powers. When arcing no longer occurred at these points, it was observed 
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to occur at randomly located places, becoming more intense at higher 
levels of average power. This observation was a factor in generating 
the hypothesis that particulate matter was inducing breakdown. The 
average power limitation imposed by the particulate matter was 
experimentally found to be lessened by allowing the arcing to continue 
until stopped, a procedure which came to be known as “processing.” 

As shown in Fig. 1, when a high average power is present, very small 
particles can cause arcing at levels well below the inherent peak power 
handling capability of the system. Processing is accomplished by apply- 
ing power to the assembled system, increasing the power levels until 
breakdown is induced by the particles and then waiting until arcing 
subsides. Processing arcs change the physical configuration of the 
initiating particles by (7) breaking them into smaller particles, (72) 
removing them from the high field region via breakdown pressure 
waves, and/or (777) vaporization. Figure 1 shows that, as the particles 


WAVEGUIDE BREAKDOWN EFFECTS 2073 


become smaller, a higher peak power can be handled before arcing 
occurs. Therefore, by raising the power level in discrete steps and 
remaining at each level for a period of time sufficient to allow processing 
to occur, the power-handling capability of the microwave system in- 
creases to a higher value. This value can be maintained as long as the 
system is intact, that is, not opened. Disassembly always introduces 
particulate matter and some reprocessing is required to achieve the 
former stable power handling level. 

Figure 5 shows an experimental processing sequence. The data was 
taken in a 10-foot ring set up as shown in Fig. 2. Parameters for this 
experiment were: 

Pulse length—120 microseconds 

Dielectric gas—SF, at 25 psig 

Ambient temperature—28°C. 

The charts appearing in Fig. 5 depict the number of ares which 
occurred in two-minute intervals at the indicated power levels. In all 
tests, duty cycle was varied by holding the pulse length fixed and 
changing the pulse repetition rate. Arcing began during Run | at a 
peak power level of 12.6 megawatts with a 1.5 percent duty cycle. It is 
noted that (7) each time the power levels were raised to higher values, 
breakdown occurred, and (72) at each power level, the number of arcs 
which occurred per two-minute interval decreased with time. At and 
above the peak power level of 66.7 megawatts, the average power was 
held constant at 1 megawatt. Performance became unstable at the peak 
power level of 139 megawatts. 

Run 2 was then executed. It is noted that the processing arcs which 
occurred during Run 1 improved system performance to the extent that 
during Run 2, no discharges occurred below the power levels of 66.7- 
megawatts peak, 1-megawatt average. Also, fewer arcs occurred at the 
higher power levels than during Run 1. Performance became unstable 
at a power level of 151 megawatts peak. The results of this experiment 
demonstrate that processing a system to a high power level raises the 
power level at which subsequent discharges occur and reduces their 
severity. 

Ideally, processing should be performed at power levels exceeding 
system requirements so that operation at the design power levels would 
be virtually arc-free. Such a procedure cannot be applied to a typical 
system which normally operates at the maximum available power level. 
Processing can still be applied in such a situation, however, if it is 
performed with the dielectric strength of the gaseous fill somehow 
reduced. This can be accomplished either by introducing a low dielectric 
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Fig. 5—Waveguide processing. 
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strength gas before processing, or by processing with the pressure of 
the normal dielectric gas reduced. Processing can then occur to a greater 
extent with the available power and will result in stable performance 
when the system is operated with normal dielectric strength. 

Dry nitrogen was first investigated as a processing agent because of 
its low dielectric strength. It was found that after the initiation of a 
processing arc, nitrogen does not recover its dielectric strength sufh- 
ciently fast to prevent breakdown on all successive pulses. Hence, 
processing cannot be satisfactorily performed using this gas. 

Processing with low pressure SF’, was next investigated. A 32-foot 
resonant ring was selected for this purpose. Results are shown in Fig. 6. 
The ring was processed to power levels of 20 megawatts peak, 700 
kilowatts average while filled with SF, at a pressure of 5 psig. After 
increasing the gas pressure to 25 psig, the peak power level was raised 
with the average power held constant; breakdown did not again occur 
until a peak power level of 35 megawatts had been achieved. Hence, 
processing with 5 psig of SF; provided almost 3 dB of margin in power 
handling capability when the system was repressurized to 25 psig. 
This data was successfully repeated after disassembly and reassembly 
of the structure. Hence, reduced pressure SF, can be employed in the 
processing of a practical system to assure some power handling margin. 


VI. PROCESSING CONSIDERATIONS AND LIMITATIONS 


The energy dissipated by processing arcs must be kept sufficiently 
low to prevent the temperature rise of the waveguide surface in the 
vicinity of the are from exceeding the melting point of the waveguide 
material. If melting of the waveguide wall occurs, the turbulence of 
the arc can eject small masses of molten metal into the waveguide 
interior, preventing processing since the additional particles induce 
continuing breakdown. Hence, the arc rate, instead of subsiding with 
time, increases. Such a situation is referred to as massive breakdown 
and is evident in Fig. 5, Run 2 at the 151-megawatt level. 

It will be shown that the power level at which massive breakdown 
occurs is inversely proportional to the square root of the time duration 
of a single arc. To analyze the effects of are duration upon waveguide 
processing, it is necessary to determine (z) the formative time of a 
processing arc, and (zz) the percentage of incident power absorbed by 
a growing arc as a function of time. The processes contributing to the 
growth of an arc which begins at the surface of a small particle have 
not, to the authors’ knowledge, been investigated in any depth. In our 
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Fig. 6—Low-pressure processing. 


investigation, no rigorous analytical determination of the desired two 
quantities was attempted; rather, both were determined experimentally. 
The experimental procedure is contained in Appendix B. The waveforms 
observed suggest that the power absorbed by a processing arc when the 
incident power is constant can be approximated by the form: 


WAVEGUIDE BREAKDOWN EFFECTS 2077 


P, = ofl — eR. (5) 
where 


P, = incident power 
T = formative time constant of the arc, independent of power level 
a = steady-state absorption percentage. 


For WR-284 waveguide pressurized to 25 psig with SI, , it was found 
that 7 = 2 microseconds and a = 1.5 percent. 

The rise in temperature of the waveguide wall in the vicinity of a 
processing discharge can be determined by assuming that the thermal 
conductivities and heat capacities of the dielectric gas and the waveguide 
material are such that power dissipated within an arc is transferred 
instantaneously to the waveguide walls. Consider the semi-infinite 
block of metal shown in Fig. 7. Let P4(¢) be the power per unit area at 
the metal surface z = 0. It can be shown’ that the temperature rise 
within the metal for ¢ > 0 is given by: 


AT = x P4(u) eis] du (6) 


2 Jo Vi-—wu 


where u is a dummy variable of integration and y is a constant dependent 
upon the thermal conductivity and heat capacity of copper. At the 
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Fig. 7—Model employed to calculate AT. 
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surface of the metal, x = 0, this integral becomes: 
=) = ui Se = du, (7) 


The power absorbed by an arc is given by eq. (5). Consequently, by 
substituting (5) in (7), the temperature rise of the waveguide surface 
at time ¢ after the arc is initiated is found to be: 


afl — e'’")P,/A [A 4, 
~ 2 a tau —uU ®) 


where A is the cross-sectional area of the arc. 


For our situation, the arc duration, t, exceeds the arc formative time 
constant, 7', so that: 





eee eee 
AT(o, t) %& a a du (9a) 
or 
AT = 7 P, Vt. (9b) 


Therefore, the temperature rise experienced at the waveguide wall 
varies as the square root of are duration, t. To prevent the waveguide 
wall temperature from exceeding the melting point of the waveguide 
material, the arc duration, ¢, must be kept small. 

To test the validity of the above theory of processing limitation 
resulting from particle generation, processing was performed within 
WR-284 resonant structures of various lengths, pressurized with SF, 
at 25 psig. The duration of an are which occurs within a resonant struc- 
ture is determined by the time required to dissipate all of the energy 
which was stored within the resonator prior to breakdown. Since stored 
energy increases with resonator length, arc duration also increases 
with resonator length. 

Whereas the incident power in a terminated line is a constant, the 
. power incident upon an arc in a resonant structure decreases with time 
as the stored energy is consumed by the arc. An expression for this 
incident power as a function of time P,(¢) is derived in Appendix B 
(eq. 60). The power absorbed by the arc is therefore: 


P,(t) = aPp(1 — e7'’") exp ite [te - TA — my} (10) 
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where 
Py, = resonator equivalent power level prior to breakdown 
v, = group velocity of electromagnetic wave within the resonator 
1 = resonator length. 


Substituting this value of absorbed power into equation (6) yields: 


ee i a(1 — e”’*) exp {—(v,/lolu — TU — &”"))} Pr, (11) 
2A Jo Vi-u 

Experimental verification of these relationships was accomplished 
on four resonant structures. The first two were waveguide cavities with 
input and output irises designed for optimum power multiplication. 
The other two were resonant ring structures. All structures were 
pressurized with SF, at 25 psig and processed as described above to 
the point of massive breakdown. The results are as shown in Table I. 

Resonators 1 and 3 were operated at a constant 0.5-percent duty 
cycle. Resonators 2 and 4 were operated with the average power held 
constant at 1 megawatt. The massive breakdown power reached is noted. 
The last column is the calculated value of the relative wall temperature 
reached at massive breakdown using eq. 11. The calculations were 
performed numerically on a digital computer and are plotted in Fig. 8. 

These calculated maximum values are relatively close, in spite of 
the fact that the structures are radically different in dimensions. The 
calculated wall temperature rise of the short resonator deviates most. 
This is probably because the power level at which breakdown occurred 
was relatively close to the theoretical breakdown limit of the waveguide 
which is 900 megawatts. These experimental data imply that massive 
breakdown occurs when the rise in waveguide surface temperature 
above ambient operating temperature exceeds a critical value of approx- 
imately 850 y/A. 

A microscopic examination of some of the parts used in these experi- 
ments confirmed the analysis. Arc marks were found to be present upon 





TABLE I—BREAKDOWN Power LEVEL VERSUS LENGTH 


Massive Breakdown 


Resonator Length Power Level AT 
No. (Meters) (Megawatts) Maximum 
] 0.135 650 675 y/A 
2 3.05 150 875 y/A 
3 3.66 150 895 7/A 
4 9.65 70 960 7/A 
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Fig. 8—Waveguide surface temperature rise during massive breakdown arcs. 


the waveguide flange gaskets which were removed from the 10-foot 
and 38-foot resonant rings after massive breakdown. One such arc mark, 
magnified approximately 200 times is shown in Fig. 9. To achieve proper 
lighting, it was necessary to tilt the gasket under the microscope. 
Hence, only a thin horizontal strip is in focus. It is noted that the arc 
mark consists of an area of discoloration approximately 6 millimeters 
in diameter, and that a small hemispherical crater approximately 
30 microns in diameter is located near the center of the discoloration. 
The area of discoloration is believed to have been caused by excessive 
heating, and the crater formed after the waveguide wall had melted. 
Referring to Fig. 1, it is obvious that the size of the conducting particle 
removed from the crater is capable of initiating additional discharges. 
The size of the particle generated is consistent with values reported 
for particles formed by de discharges of roughly the same energy as 
experienced at massive breakdown in the resonators.° 

Processing considerations in a practical system based upon the above 
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theory of particle generation and massive breakdown can be illustrated 
by the following example. Suppose a microwave system is to be designed 
to transmit 20 megawatts of peak power with 60-microsecond pulses, 
and that after 10 microseconds of a pulse has elapsed, an arc is initiated 
by a particle of suitable size. The arc duration is therefore 50 micro- 
seconds, and eq. (9b) predicts a temperature rise of: 


AT = 7 X 20 X 10°V50 X 10° = 141 x 10°F 
For SF, , a = 1.5 X 10°? => AT = 2100 y/A. Hence, melting occurs, 
and massive breakdown will ensue. If, however, the initiation of the 
are was detected (such as by comparing forward power levels at the 
amplifier and at the load) and the RF drive was blanked within, say, 
4 microseconds after arc initiation, then: 


A 


AT = % x 20 x 10°W/4 X 10°° = 40 X 10°ay/A. 





Fig. 9—Are mark. 
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Fora = 1.5 X 10°’, AT = 6007/A. Thus, by employing this “early-off” 
protection, the critical temperature is not exceeded, and processing can 
be carried to completion in a high average power system. 


VII. CONCLUSION 


The theoretical and experimental investigation of the effects of high 
average power levels and long pulse lengths on waveguide breakdown 
has shown that: 

(z) Loose particulate matter (e.g. thermally detached from the 
waveguide wall) is heated to high temperatures by high levels 
of average power. The resulting rarification of the dielectric gas 
in the vicinity of the heated particle, coupled with field enhance- 
ment at the particle boundary, causes breakdown to occur at 
values of peak power below that expected in the absence of 
particles. 

(it) Waveguide “processing”? with RF discharges of short duration 
can be used to eliminate undesirable foreign matter and increase 
the power handling capability of a microwave system. 

(iii) To achieve processing, arcs must terminate as soon as possible 
after they occur by blanking the remainder of the RF pulse. 
This prevents excessive energy dissipation in the arc which can 
cause waveguide surface melting. 

(tv) The ultimate peak power level to which a system can be processed 
varies inversely with the square root of the allowed arc duration. 

(v) Processing requires a dielectric gas which rapidly recovers its 
dielectric strength. Sulphur hexafluoride is ideal and perhaps 
the only practical choice. With its use, processing at reduced 
pressure can achieve some breakdown margin in the usual case 
where processing power is limited to the normal operating power. 
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APPENDIX A 


Field Enhancement and Heating of Conducting Particles 


Consider a spherical conducting particle located within a waveguide 
as shown in Fig. 10a. For simplicity, the particle is assumed to be 
centrally located and completely detached from any of the waveguide 
walls. When electromagnetic waves are propagated through the wave- 
guide, the boundary conditions imposed by the particle cause the field 
to be distorted. Since the radius of the particle is much less than a 
wavelength, the quasi-static approximation can be applied to obtain 
a good estimate of the fields existing in the vicinity of the particle. 
The models shown in Figs. 10b and 10c will be used to obtain the quasi- 
static fields. The coordinate systems used in Figs. 10b and c are un- 
related to each other or to that of Fig. 10a. The zero-order longitudinal 
magnetic field can be neglected, since the particle is assumed to be 
centrally located (H, is zero at the center of the waveguide transverse 
section for the mode of interest, namely the TE,)). 





(b) (c) 


Fig. 10—Particle model: (a) Particle location. (b) EF field model. (c) H field model. 
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The quasi-static approximation consists of expressing the total 
electromagnetic field in the form: 


E=e +e, +we,+-::- (12) 
H =h, + oh, + wh, + -:: (13) 
where 
VxXe=0 (14) 
V Xho =0 (15) 


oY X hy = +e 5 


oh 
woV x e; = —Ho eva 

0Cn-1 
woV Xh, = +60 ag (16) 
wT X @, = — ho Te (17) 


The approximation consists of truncating (12) and (13) after a sufficient 
number of terms. The zeroth order fields ep and hy are found from eqs. 
(14), (15) and the unperturbed fields EZ, and Hy. 

The magnitudes of the higher-order correction fields are related to 
the magnitudes of the zeroth order fields by the factor (a/\)”, where 
\ is the wavelength, a is a characteristic dimension of the perturbation, 
and 7 is the index of the correction field being considered. Since a < )}, 
it follows that the first and all higher-order correction fields can -be 
neglected, and the fields in the vicinity of the particle are closely 
approximated by the zeroth-order terms. These are found by solving 
LaPlace’s equation in spherical coordinates for the electric and magnetic 
scalar potentials, taking the gradients of these functions, and applying 
the boundary conditions € |, = —oiz 3 €oy |p-a = 0, Bo |p = Aoi: ; 





ho, |-=a = O. The results are: 
a\* a\* 
e&) = -B,|1 + 2(*) | cos #i, + Bla ~- (2) | sin Big (18) 
a\* a\* 
hy = il 1 — (2) | cos 6i, — 1] 2 + (2) | sin Gig . (19) 


WAVEGUIDE BREAKDOWN EFFECTS 2085 


Now, from eq. (18), it follows that e, assumes its maximum value 
at 6 = 0,r = a. This value is: 


| Co [ape a 3k, * (20) 
Hence, the field enhancement factor, 8, for spherical conducting particles 
is 3. 
From eq. (19), 

hy |-2 = —3H, sin Big (21) 

| I, | a | hy(a) | = 3H, sin 6 (22) 

where J, is the surface current density induced in the particle by the 

field. If the microwave frequency is high enough such that the skin 

depth of penetration into the particle is much less than the particle’s 

radius, the surface current density may be assumed to flow uniformly 
within one skin depth, 6. Hence, the current density J becomes: 


I, 
ee (23) 


The power dissipated per unit volume within one skin depth of the 
particle is therefore: 


Be sedge ae (24) 
~ 2 ¢ 2 68 


where o is the conductivity of the particle. The power dissipated per 
unit surface area of the particle therefore becomes: 


eee eal (25) 


2 206 
But 
6 = V 2/wuoo (26) 
p, = LET, fave 27) 
2 20 
or 


_ 9Ho sin’ 4 feopto 
a 2 Qo 28) 


The total dissipated power becomes: 


Pun = | Px dS (29) 


wT 2a 2 
= i i 9H Ae a’ sin® 6 dé dé (30) 
0 0 2 20 
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or 
_ 36H2a° ee 
ee er 2 Qe : (31) 
But the propagating power level, P, is related to H, by the expression 
P = ae ? (32) 


where Z, is the TE,, wave impedance. 


14470’. alae 
ed a 


Pains = 3Zd,d> 


(33) 
Since 


Uo w/c 


& Ve! — (n/a) ' 
14470” | 2 
P iss — at _ (<=) lp 35 
3did,. NQo E dyw (85) 
where c is the free-space velocity of light. 


To calculate the temperature rise of the particle resulting from this 
dissipated power, let n be the heat flux density flowing from the particle: 


= 


(34) 





$ MR Pac (36) 


particle 
surface 


Now, assuming heat flow by conduction through the gas, 
n= —§EVT (37) 


where £ is the thermal conductivity of the gas. In the steady-state, 
conservation of energy implies that: 


Vn =0 (38) 
—V-[EVT] = 0. (39) 

For SF, ,” 
2,7 (40) 


where £, = 4.3 X 107° watts/meter-°K’. Hence, 
EV -[TVT] = 0. (41) 
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From circular symmetry, eq. (41) becomes 


bd ap aE 
2) ord = 0. (42), 


The solution of this differential equation is: 


| Ci 
= Cg — Er (43) 


where ¢;, C2 are constants determined by applying the boundary 
conditions 


Pim =T., | dA = Pa. 


particle 
surface 


= 2 Praise. 
T= Ti+ 53 (44) 


The absolute temperature of the particle, T, , is therefore 


The result is: 


Paves 


— = 2 
di nes T ee dies =e Qrkoa 


(45) 
where P,;,, is given by eq. (85). It is observed that P4;,, is dependent 
upon the electrical conductivity, o, of the particle. 

Since the electrical conductivity of a conductor varies inversely 
with the absolute temperature, 


rT, 
= 6; T, (46) 
where 7’, is the ambient temperature and is assumed to be 293°] and 
o, 1s the conductivity at ambient temperature. For copper, the bulk 
conductivity is equal to 5.5 X 10’ meters/ohm. Because microwave 
current densities flow within a narrow skin depth, ripples in the surface 
of the conductor cause a decrease in the effective conductivity near the 
surface.” At S-band frequencies, surfaces with rms variations on the 
order of several skin depths (such as rough spheres) have an effective 
conductivity of approximately 62 percent of the bulk conductivity. 
Hence, a value of 3.44 X 10’ meters/ohm will be assumed for o, when 
calculating the temperatures of the copper spheres. 


2088 THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1972 


Substituting eqs. (35) and (46) into (45) yields: 


=e) 
2 d yj 


T= 2 
p T. —— : 47 
= 3d,d,£, Vo,T 4 ( ) 


724 | (1 7 =! 
2 dw 


B = ——_ 48 
3d, do, V Ool'4 ( ) 


Letting 





yields: 
T? = T? + BaP, VT, . (49) 
Experimentally, d; = 2.84 inches, d, = 1.34 inches. 
B = 4.68 X 10’ in MKS units 


Finally, 7, is equal to the temperature of the waveguide in the region 
where the particle is located. Experimentally, the particles were intro- 
duced at the waveguide flanges, and the flanges were found to experience 
a temperature rise of 65°C/megawatt of average power above ambient. 
Hence, to be consistent with experiment, 


T= 7, +65 X10 “P., (50) 
where 7', is the ambient room temperature (293°K). Hence, 
T; = (293 + 65 X 10°°P,,)” + 4.68 X 10°aP.. VT). (51) 


Values of 7’, which satisfy (51) appear in Table II for various values 
of P,, and particle radius a. 


APPENDIX B 


Arc Formative Time and Power Absorption 


The formative time of a processing arc is defined as the interval from 
arc initiation to complete interruption of power transmitted. Photo- 
graphs of traveling-wave resonator waveshapes taken during the 
occurrence of processing arcs indicate that processing-arc formative 
times are independent of peak power level over the range of investiga- 
tion (11 to 80 megawatts peak). 

The formative time was measured within a terminated WR-284 
system pressurized with SF, at 25 psig. Transmitted power pulses were 
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TaBLe I]—ParticLe TEMPERATURE T7', AS A 
FUNCTION OF S1zE AND PowrErR LEVEL 


Average Power (Watts) 
Particle Radius | |————-—————, — J 


(Microns) 105 2x 105 | 410° | 8 x 105 | 1.6 X 108 
60 388°K | 478°K | 648°K | 960°K 1425°K 
30 346°K | 395°K | 491°K | 664°K 964°K 
15 323°K | 350°K | 405°K | 508°K 696°K 

7.5 310°K | 327°K | 365°K | 426°K 550°K 
2 300°K | 312°K | 332°K | 368°K 436°K 


monitored on a storage oscilloscope at the termination. The formative 
time was found from the rate of decay of the transmitted power pulse 
whenever a discharge occurred at any point between the final power 
amplifier and the termination. A terminated (non-resonating) structure 
was employed for this measurement to divorce the arc formative time 
from the natural decay of a resonator caused by the additional arc 
loading. It was observed that the rate of decay of power transmitted 
past an arc was independent of the initiating power level and that the 
transmitter power waveshape decayed exponentially with a time con- 
stant of 2 microseconds. Hence, in WR-284 waveguide pressurized to 
25 psig with SF, , the power absorbed by an arc can be written as: 


Pi = aP,[1 = et) (52) 
where 
t = time 
P, = incident power 
a = steady-state percentage of incident power absorbed by arc 
T = 2 microseconds. 


The value of a was determined by use of a traveling wave resonator. 

When an arc is struck within a resonant structure, it begins to absorb 
the energy that was stored in the resonator and the fields decay to zero. 
The rate of decay of the fields is dependent upon the percentage of 
incident power absorbed by the arc. This relationship is derived as 
follows. 

Let Pr(t) be the instantaneous power level within the resonator; 
then, when an arc is initiated, 


Prt + 7) = [1 — a@]Pr@) (53) 
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where 


I 


transit time through the resonator (= l/v, where I is the 
resonator length and v, the group velocity) 

a(t) = percentage of incident power absorbed by arc as a function 
of time. 


T 


Since 7 is small, 
dP, 





Prt + 7) & P(t) + 7 Ts (54) 
Substituting eqs. (52) and (54) into (53) yields 
‘- o es bai op: (55) 
=P, = C exp {—e/7[t + Te’ J}. (56) 
Now, 
P,(t = 0) = Pr, (57) 
=> P, = Pp, exp {—a/7[t — TIL —e’")]}. (58) 
But 
Rage gee a 59) 
fo aes ts OTe | ( 
Py, & Pp, exp {-—a/7[t — TA —14+ t/T — #/2T)}} (60) 
= P,, exp {—a/7[t — t+ #/2T]} (61) 
= P,(t) & Pr, exp {—at’/2Tr}. (62) 


Figure 11 shows the typical decay of resonator forward power during 
the occurrence of an arc. By fitting eq. (62) to this curve, the value 
of a can be found since the ring length 1, the group velocity v, , and 
the formative time constant 7’ are known. The value of a was experi- 
mentally found to be 


w= 0015, | (63) 


Hence, in the steady-state, an are struck in WR-284 waveguide filled 
with 25 psig of SF; absorbs 1.5 percent of the power which is incident 
upon it and reflects the remainder. 

It is noted that the above method for determining the percentage of 
power absorbed by an arc provides greater accuracy than could be 
obtained by direct measurement of the power levels incident upon and 
reflected from the arc since the latter involves the subtraction of two 
nearly equal quantities. 
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10 FOOT RESONANT RING 
DIELECTRIC—SF—25 P.S.I.G. 


PEAK POWER IN MEGAWATTS 





0 0.5 1.0 1.5 2.0 2.5 3.0 
TIME IN MICROSECONDS 


Fig. 11—Resonant ring power decay. 
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A Table Look Up Approach to Loop Switching 


By L. H. BRANDENBURG and B. GOPINATH 
(Manuscript received July 24, 1972) 


In this paper we consider some questions of implementation of a 
scheme described in Section V of Ref. 1 for addressing message blocks 
in the Pierce loop system.” The scheme consists of using a stream of 
binary digits (0 and 1) as the address of the destination loop of a message 
such that the scalar product* of the address with a stream of binary 
digits of equal length stored at a loop gives the distance between the 
loop and the destination. The message is routed along a path that 
minimizes distance between source and destination. The binary streams 
used in this scheme can be obtained by factoring the distance matrix D 
of the graph representing the connection of loops into two binary- 
valued matrices P and Q such that 


D = PQ, 


where the superscript ‘‘?’’ stands for matrix transposition. The kth 
row of P represents the address to be prefixed to a message destined 
for loop k. The kth row of Q represents a binary sequence to be stored 
in loop k. The scheme discussed in, Section V of Ref. 1 provides a par- 
ticular way of factoring D. For completeness, we give a description of 
that factorization as follows. 

Let s be the diameter of a graph on n vertices and let D be its distance 
matrix. (Note: s = max,; D;;). We consider a PQ‘ factorization of D 
with P and Q matrices each of order n X ns. The 7th row of P, P; , is 
zero everywhere except in positions (s(¢ — 1) + 1) to sz where it is one. 
The ith row of Q, Q; , is constructed as follows: for every j, there are 
exactly D;; 1’s in positions (s(j — 1) + 1) to sj with the rest of the 
entries of Q; equal to zero. 


* The scalar product of two binary streams of equal length is the number of places 
they both have a “‘one’’. 
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s 2s 3s ns 
J vy 
1---1 0---0 0---0 
pe, 0:---0 1---1 0 
0 0 0:---0 O 0 1 1 
8 Diz 2s Drs 38 Ds. ns 
1 i 1 i 
Oh SPeser ps ka Wee a eee. {iets Ee eee 6 
1---10 0 1---1---0 1 1 0 
Q= 
fete le ated, ae Leese Tae 0 re 0 
Dar Dus D3 


Obviously the length of addresses in this scheme is ns. The P matrix 
defined above, rows of which are addresses to be prefixed onto messages, 
is the same for all graphs with diameter s. The Q matrix, rows of which 
are stored in the loops in the network, contains the information that 
identifies a particular graph. 

As will be discussed further, a simple mechanization of this scheme 
reveals that it is essentially an obvious table look up scheme. In a table 
look up scheme, each loop has stored the distance between it and every 
other loop in the networks, and each address essentially provides a 
signal to look up or read out a particular distance that is stored. 

The above PQ‘ factorization can be transformed into a table look 
up scheme in the following way. The integer 7 (Sn) uniquely specifies 
the 7th row of P; thus, a message address need only consist of the 
binary representation of 7, which, of course, requires at most [log. n]* 
bits. Each loop, instead of storing simply a row of matrix Q, stores 
a binary array consisting of n rows and [log. s] columns. The 7th row 
of this array contains the binary representation of the distance between 
the loop in question and the ith loop. The binary array can be im- 
plemented by a read only memory the 7th row of which is accessed by 
the binary sequence representing integer 7. Note that in the original 
PQ’ form of this scheme, both matrices P and Q were dependent on s, 
the diameter of the graph. In the table look up mechanization, s appears 


* [x] means ‘“‘the least integer greater than 2’’. 
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only in the dimension of the stored memory. (The authors had originally 
considered a table look up mechanization for minimum distance routing 
but had concentrated their early efforts on the more mathematically 
fruitful schemes discussed in Refs. 1 and 3. It is interesting to note at 
this point that the table look up mechanization has evolved in a natural 
way from a special case of one of those schemes.) 

The table look up scheme meets the following important criteria: 

(2) The simplicity of constructing addresses is important in the 
case of large graphs with no special structure. Even if one were to find 
minimum length addresses similar to the ones described in Sections I, 
II, and III of Ref. 1, we suspect that the construction of these addresses 
would be very complicated. In the present case, each address is obtained 
trivially as the binary representation of an integer. 

(iz) Since present plans for length of message blocks envisage lengths 
of perhaps a few thousand bits, the length of the message address is 
an important parameter in any large loop system. The present method 
guarantees minimum length addresses, [log, J], provided that the 
graph is not constrained to have a particular structure. 

(iat) The scheme can be adapted by prefixing some control digits 
to do alternate routing. 

(wv) It is simple to determine the necessary changes in the stored 
memory required by updating (adding loops) or modifying the network, 
since such changes can be identified by inspection of the distance 
matrix. 

(v) The speed of operation of this scheme can be very fast compared 
to the speed of the message. 

(vt) If the graph is hierarchical in the sense of Pierce,” except that 
interconnections among loops at the local level are allowed, then routing 
can be accomplished by using the Pierce scheme at higher levels, and 
the present scheme at the local level. 

At first it might appear that the present scheme requires longer 
address lengths than the Pierce addressing scheme” for the case of 
strictly hierarchical graphs. For most cases, we show this is not so. 
A hierarchical graph of k levels can be represented as in Fig. 1. 

N; represents the maximum of the number of branches connecting 
each vertex at the (¢ — 1)st level, to vertices at the ith level. Pierce’s 
scheme obviously requires addresses of length > /-, [log N;] bits. The 
present scheme requires at most 


m = flog (1 + N, + NiN. + «+: + N,N, «++ N,)] bits. 
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Ist LEVEL 


2nd LEVEL 





kth LEVEL Ne 


Fig. 1—A hierarchical graph of k levels. 


But 











1 1. 1 
a [10g Nv - mi a Ni * NiNy-1 ee N.N,, a x) |: 
If, at each level, there is at least one vertex with two or more branches 
descending to the next level (this will always be the case if the levels 
“fan out’’) then N; = 2 for each 7. Then 


1 1 1 
ms |bgviv,--m(1+4+h+--- +h) 
< [log 2N,N. --- N,] S [log NiN2 --: Nz] +1. 


Therefore, if each N; is a power of 2, then Pierce’s scheme gives ad- 
dresses at most one bit shorter than those of the present scheme. How- 
ever, if even one of the N,’s is not a power of 2, then the present method 
has address lengths at most equal to those of Pierce’s scheme. In fact, 
even if each N; is a power of 2, the Pierce scheme has the one bit ad- 
vantage if and only if every vertex at the zth level has exactly N,,, 
branches connected to the (¢ + 1)st level, for all levels. 

In order to compare the present scheme with other schemes described 
in Ref. 1, we assume that the important parameters are of the following 
orders: 


s = diameter of the graph + 6 
nm = number of loops ~ 10° 
length of message block ~ 4 X 10’ bits. 
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Then for the present scheme, the address length ~ 10 bits, or 1/4 
percent of the message block length. The memory requirement at each 
loop is n log, s or ~ 3 X 10° bits. For these numbers, the presently 
available memories, either core or the more advanced semiconductor 
memories, could be used to realize very high speed routing of messages. 

For the other schemes in Refs. 1 and 3, the address length was con- 
jectured not to exceed 2(n — 1) bits = 2 X 10° bits, or 50 percent of 
the message block length. The memory requirement is of the same order 
or about 2/3 of that of the present scheme. 

With regard to time of computing distance at each junction, the 
table look up scheme using semiconductor memory might require 
on the order of 100 nanoseconds; on the other hand, the schemes dis- 
cussed in Refs. 1 and 3, based on inner product or Hamming distance 
calculations, might require on the order of 5,000 nanoseconds. 

These numbers, together with the six points listed previously, provide 
evidence of the practical advantage of the table look up scheme over 
the other schemes in Ref. 1. We have, of course, omitted discussion 
of such important questions as reliability and the implementation of 
updating the stored memory in each loop due to changes in the net- 
work. However, these questions are of equal importance in the present 
scheme and the schemes in Ref. 1, since the basic routing rule is the 
same (minimum distance) and the sizes of the stored memories are 
comparable as shown above. 

The authors are indebted to H. 8. McDonald for helpful discussions 
concerning some of the practicalities of implementing a loop switching 
system. 
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